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ABSTRACT 

We present an extension of traphic, the method for radiative transfer of ionising ra- 
diation in smoothed particle hydrodynamics simulations that we introduced in Pawlik 
& Schaye (2008). The new version keeps all advantages of the original implementa- 
tion: photons are transported at the speed of light, in a photon-conserving manner, 
directly on the spatially adaptive, unstructured grid traced out by the particles, in 
a computation time that is independent of the number of radiation sources, and in 
parallel on distributed memory machines. We extend the method to include multiple 
frequencies, both hydrogen and helium, and to model the coupled evolution of the 
temperature and ionisation balance. We test our methods by performing a set of sim- 
ulations of increasing complexity and including a small cosmological reionisation run. 
The results are in excellent agreement with exact solutions, where available, and also 
with results obtained with other codes if we make similar assumptions and account for 
differences in the atomic rates used. We use the new implementation to illustrate the 
differences between simulations that compute photoheating in the grey approximation 
and those that use multiple frequency bins. We show that close to ionising sources the 
grey approximation asymptotes to the multi-frequency result if photoheating rates are 
computed in the optically thin limit, but that the grey approximation breaks down 
everywhere if, as is often done, the optically thick limit is assumed. 
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- diffuse radiation - cosmology: large-scale structure of Universe 



1 INTRODUCTION 

New telescopes such as Planck\ LOFAR^ MWA^, ALMA^ 
and JWST^ will soon open up new windows onto the epoch 
of reionisation (e.g., Barkana & Loeb 2001; Ciardi & Ferrara 
2005; Fan, Carilli, & Keating 2006; Furlanetto, Oh, & Briggs 
2006 for reviews of this epoch). Data collected by these tele- 
scopes is expected to shed light on many unresolved issues in 
our current understanding of how galaxies form and evolve 
and interact with their surroundings. Detailed theoretical 
studies, however, will be needed to interpret it. Amongst 
the most promising techniques to perform such studies are 
cosmological simulations of reionisation. 

Modern simulations of reionisation aim to combine the 
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first-principle modelling of the gravitational growth of den- 
sity fluctuations and of the hydrodynamical evolution of the 
cosmic gas in the expanding Universe with recipes for star 
formation and associated feedback and to follow also the 
propagation of ionising radiation emitted by the first ion- 
ising sources. The computationally efficient, but accurate 
implementation of the radiative transfer (RT) is currently 
one of the biggest challenges for simulating reionisation. 

Computing the ionising intensity throughout the sim- 
ulation box requires solving the seven-dimensional (three 
space coordinates, two directional coordinates, frequency 
and time) RT equation. This is a formidable task, not only 
because of the high dimensionality of the problem, but also 
because of the large number of ionising sources contained 
in typical cosmological volumes. To accomplish it, existing 
approaches (e.g., Abel, Norman, & Madau 1999; Gnedin 
& Abel 2001; Ciardi et al. 2001; Nakamoto, Umemura, & 
Susa 2001; Maselli, Ferrara, & Ciardi 2003; Razoumov & 
Cardall 2005; Mellema et al. 2006; Susa 2006; Ritzerveld & 
Icke 2006; McQuinn et al. 2007; Semelin, Combes, k Back 
2007; Trac & Cen 2007; Pawlik & Schaye 2008; Aubert & 
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Teyssier 2008; Altay, Croft, & Pelupessy 2008; Petkova & 
Springel 2009; Finlator, Ozel, & Dave 2009; Gritschneder 
et al. 2009; Paardekooper, Kruip, & Icke 2010; Hasegawa 
& Umemura 2010; Cantalupo & Porciani 2010; Parti et al. 
2010) to transport ionising photons must often resort to a 
number of approximations. 

The accuracy of several ionising (cosmological) RT 
codes has been assessed in test simulations that were per- 
formed as part of a series of comparison projects (Iliev et 
al. 2006a; Iliev et al. 2009). The results of the comparisons 
are encouraging and indicate that the participating codes 
have reached a certain level of maturity (Iliev et al. 2009). 
The design of most of the test simulations was kept sim- 
ple in order to facilitate comparisons between different RT 
codes. More recently, the performance of different RT codes 
has been compared in cosmological simulations of reionisa- 
tion with an equally promising degree of agreement (Zahn 
et al. 2010). However, the inclusion of RT in state-of-the-art 
simulations of structure formation remains a tough compu- 
tational challenge, as we now explain. 

RT codes that are both spatially adaptive and parallel 
on distributed memory are still rare (see, e.g.. Table 1 in 
Ihev et al. 2006a and Table 1 in Ihev et al. 2009). Nearly ah 
reionisation simulations are therefore performed on uniform 
grids. Combined with the fact that large simulation boxes 
are needed to model representative volumes of the Universe, 
this means that the spatial resolution of state-of-the-art RT 
simulations of reionisation is typically far below that of the 
underlying spatially adaptive hydrodynamical simulations. 
In fact, many RT simulations of reionisation ignore hydro- 
dynamical effects altogether and assume the gas traces the 
dark matter. Small-scale structure in the cosmic gas is there- 
fore often ignored or included only in a statistically sense. 

Cosmological simulations of reionisation typically con- 
tain millions of star particles (e.g., Iliev et al. 2006b). Large 
numbers of ionising sources pose a challenge to simulations 
of reionisation because for most of the existing RT meth- 
ods the computation times increases linearly with the source 
number. The usual practice of reducing the number of ionis- 
ing sources by combining sources that fall into the same cell 
of a superimposed mesh renders reionisation simulations fea- 
sible, but also reduces the spatial resolution at which the RT 
is performed. Note that the inclusion of diffuse ionising radi- 
ation emitted by recombining ions further increases the num- 
ber of ionising sources. To reduce the computational effort, 
this recombination radiation is therefore usually treated us- 
ing the on-the-spot approximation (e.g., Osterbrock 1989), 
which assumes it to be re-absorbed in the immediate vicin- 
ity of the recombining ion. However, the validity of this ap- 
proximation remains to be assessed (e.g., Ritzerveld 2005; 
Williams & Henney 2009, Hasegawa Sz Umemura 2010). 

RT simulations of reionisation are still often performed 
by post-processing pre-computed static density fields. This 
static approximation is appropriate for simulating the ini- 
tial phase of rapid growth of ionised regions or the propa- 
gation of ionisation fronts on cosmological scales (see, e.g., 
the discussion in Iliev et al. 2006b). Once the speed of ioni- 
sation fronts becomes comparable to the sound speed of the 
ionised gas, the static approximation, however, becomes in- 
applicable and a full radiation-hydrodynamical treatment is 
required. In any case, the static approximation breaks down 
after about a sound-crossing time, as the Jeans filtering of 



the gas can then no longer be ignored (e.g., Gnedin 2000). 
Although radiation-hydrodynamical feedback from reionisa- 
tion is known to play a key role, most of the large-scale 
reionisation simulations performed to date ignore it. 

In Pawlik & Schaye (2008, hereafter Paper I) we pre- 
sented the RT method TRAPHIC (TRAnsport of PHotons 
In Cones) for use in Smoothed Particle Hydrodynamics 
(SPH; Gingold & Monaghan 1977; Lucy 1977) simulations. 
TRAPHIC can be used to solve both the time-independent and 
the time-dependent RT equation in an explicitly photon- 
conserving manner. It employs the full spatial resolution of 
the underlying SPH simulation because it works directly on 
the unstructured grid formed by the discrete set of SPH par- 
ticles. It achieves directed transport of radiation on the irreg- 
ular distribution of SPH particles by tracing photon packets 
inside cones. The solid angle of these cones thereby sets the 
angular resolution at which the RT is performed. TRAPHIC is 
by construction parallel on distributed memory machines if 
the SPH simulation itself is parallel on distributed memory 
machines. 

The computational cost for simulations with traphic 
is independent of the number of ionising sources. It merely 
scales with the product of the number of spatial and angular 
resolution elements, i.e. with the number of SPH particles 
and the number of cones needed to tessellate the sky. For 
comparison, the computational cost of conventional ray and 
photon tracing methods scales with the product of the num- 
ber of spatial resolution elements (gas particles or gas cells) 
and the number of sources. Since the number of sources is 
typically proportional to the number of spatial resolution 
elements, conventional ray and photon tracing methods face 
an expensive scaling with the square of the number of spatial 
resolution elements. In contrast, a relatively small number 
of angular resolution elements is typically sufficient to ob- 
tain converged results (Paper I; see also, e.g., Trac Sz Cen 
2007; Paardekooper, Kruip, & Icke 2010), and this, in com- 
bination with the independence of the computational cost 
of the source number, makes TRAPHIC ideal for simulations 
containing large numbers of sources (as is the case for, e.g., 
reionisation simulations) as well as for an explicit treatment 
of the diffuse radiation component. 

In Paper I we presented an implementation of TRAPHIC 
for use on (sets of) static density fields in the SPH code GAD- 
GET (Springel 2005). We applied this implementation to the 
transport of monochromatic ionising radiation in hydrogen- 
only gas at a fixed temperature. We demonstrated its excel- 
lent performance in several (static density field) test prob- 
lems that were designed to allow a detailed comparison to 
results obtained with other RT codes. Here we describe, test 
and discuss an extension of this implementation. The new 
implementation of TRAPHIC allows for the transport of multi- 
frequency radiation in primordial gas, i.e., in gas consist- 
ing of both hydrogen and helium. In addition to the com- 
putation of the ionisation state, it also allows for the self- 
consistent computation of the temperature of photoionised 
gas. The new implementation still solves the RT equation 
only on static density fields. The radiation-hydrodynamical 
coupling of TRAPHIC will be described in a future work. 

Because the computational cost for RT simulations is 
typically proportional to the number of frequencies at which 
the RT equation is solved, many RT simulations of ionis- 
ing radiation discretize the RT problem using only a single 
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frequency bin. In the grey approximation (e.g., Mihalas & 
Weibel Mihalas 1984), ionising radiation within this bin is 
then assigned an effective absorption cross-section and an ef- 
fective photoenergy that is injected in the gas upon absorp- 
tion (for examples see, e.g., Iliev et al. 2006a). The corre- 
sponding photoionisation heating rates can be computed as- 
suming the optically thin or the optically thick limit. We use 
our new implementation to show that RT simulations which 
employ the grey approximation yield gas temperatures that 
agree with the exact multi-frequency solution only if grey 
photoheating rates are computed in the optically thin limit 
and then only close to the ionising sources. Grey RT sim- 
ulations that compute photoheating rates in the optically 
thick limit significantly overestimate the typical tempera- 
tures of the photo-heated gas. A related treatment of this 
subject and a discussion of its astrophysical implications can 
be found in Abel & Haehnelt (1999). 

The structure of this paper is as follows. In Sec. 2 we 
present a brief review of the main concepts behind TRAPHIC. 
In Sec. 3 we then discuss the equations that govern the evolu- 
tion of the ionisation state and temperature of gas exposed 
to ionising radiation. With these preparations in hand we 
are ready to present our new implementation of TRAPHIC in 
Sec. 4 (and in the appendix). We discuss an extensive set 
of tests of this implementation (on static density fields) in 
Sec. 5. There we also investigate the applicability of the grey 
approximation by comparing simulations using the grey ap- 
proximation with simulations using multiple frequency bins. 
We conclude with a brief summary in Sec. 6. 



2 TRANSPORT OF PHOTONS IN CONES 

In this section we briefiy summarise the basic concepts un- 
derlying the RT method TRAPHIC and introduce some of the 
notation that will be frequently employed in the following 
sections. The reader is referred to the original description 
in Sec. 4 of Paper I for more details as that description re- 
mains valid. The extensions presented in this work concern 
only the application of TRAPHIC to the transport of ionising 
photons, which will be described in Sec. 4. 

To introduce essential notation we briefiy recall that 
SPH is a Lagrangian numerical method to solve the Euler 
equations of fiuid dynamics through the representation of 
continuum fiuids by discrete sets of particles (for reviews 
see, e.g., Monaghan 2005; Springel 2010). Any property, say 
Aj, of any given particle i is determined by performing a 
weighted average, or smoothing, Ai — TUj / pjAjWij{h) 
of the corresponding property Aj of all other particles j, 
where rrij and pj are the mass and density of particle j, 
and Wij(h) is the SPH kernel that depends on the SPH 
smoothing length h. 

In the following we make the common assumption that 
the kernel Wij is compact so that there is a finite number 
A^ngb of neighbouring particles j within a sphere of radius h 
around particle i. If the smoothing lengths, which represent 
the spatial resolution elements of the SPH simulation, are 
allowed to vary in space such that the number of neighbours 
A^ngb remains fixed, then SPH enables simulations whose 
spatial resolution adapts to the fiuid geometry. It is this 
feature that makes simulations with a large dynamic range 
possible and that is perhaps the main reason for the numer- 



ous and successful applications of SPH to solve multi-scale 
astrophysical problems such as galaxy formation and reion- 
isation. 

TRAPHIC transports photons directly on the SPH parti- 
cles (i.e. without interpolation to a superimposed numerical 
grid) and hence the full dynamic range of the SPH simula- 
tion is employed. The photon transport can be decomposed 
into the emission of photon packets by source particles fol- 
lowed by their directed propagation on the irregular set of 
SPH particles. We now briefiy describe both these parts in 
turn. 

Photon packets, each of which carries photons of a char- 
acteristic frequency ly, are emitted from source particles to 
their iVngb neighbouring SPH particles (residing in a sphere 
of radius h centred on the source) using a tessellating set 
of A'^c emission cones. The number of neighbours iVngb is a 
parameter that determines the spatial resolution and is usu- 
ally matched to the number of neighbours Nngh (residing in 
the sphere of radius h) used in the computation of the SPH 
particle properties, i.e., iVngb ^ A^ngb- The number of cones 
A^c is a parameter that determines the angular resolution of 
the RT. Emission cones are necessary to achieve an isotropic 
emission despite the generally highly irregular distribution 
of SPH particles (see App. A in Paper I). The last param- 
eter that controls the emission of photon packets by source 
particles is the number Ni, of frequency bins that are used 
to discretize the associated radiation spectrum. 

Each of the emitted photon packets (of frequency ly) has 
an associated propagation direction. This propagation direc- 
tion is chosen to be parallel to the central axis of the corre- 
sponding emission cone. After emission, the photon packets 
are traced downstream along their propagation directions. 
The packets thereby remain confined to the solid angle of 
the emission cone into which they were originally emitted 
thanks to the use of transmission cones with solid angle 
A-k/Nc. The transmission cones prevent the unconfined dif- 
fusion of photon packets on the unstructured grid of SPH 
particles that would otherwise occur and they ensure that 
the transport is directed. Because the transmission cones are 
defined locally (i.e., at the positions of the transmitting par- 
ticles) and because photon packets are only transmitted to 
the subset of the TVngb nearest neighbours that are inside the 
transmission cones, the angular resolution at which photon 
packets are transferred is independent of the distance from 
the source (even though the surface areas implied by the 
solid angles of the original emission cones increase with the 
square of the distances from the corresponding sources) . As 
a result, the sharpness of the shadows cast by opaque objects 
such as dense neutral clumps and filaments is independent 
of their distances from the sources. 

Virtual particles (ViPs) are introduced to accomplish 
the photon transport along directions for which no neigh- 
bouring SPH particle in the associated emission or trans- 
mission cones could be found. The properties of the ViPs 
(like, e.g., their densities) are determined through SPH in- 
terpolation from the iVngb neighbouring SPH particles. Their 
name refers to the fact that ViPs are temporary constructs 
that are only invoked to accomplish the transport of photon 
packets in empty cones. They do not hold permanent infor- 
mation, they do not affect the SPH simulation and they are 
deleted as soon as they have fulfilled this task. 

The photon transport is supplemented with a photon 
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Figure 1. Flow charts. Left chart: overview of traphic. The RT simulation starts with the assignment of emissivities to source particles. 
Thereafter, photon packets are transported from all particles to their neighbours (see right chart for details). ViPs, if needed, are created 
in advance of the transport and are deleted immediately afterwards. The transport cycle continues until the user-defined stopping criterion 
is satisfied. Then, the properties of the particles are updated according to the photon-gas interactions that occurred. Finally, the RT time 
step is advanced. The RT simulation ends at time Tr. Right chart: details of the transport of photons to neighbours (grey box in the left 
chart). Photons in photon packets are distribute amongst the iVngb neighbouring particles inside A^c (tessellating) emission or (regular) 
transmission cones with solid angle Att/Nc centred on the corresponding propagation directions. The optical depth to each neighbour is 
computed and the fraction of transmitted and absorbed/scattered photons is determined. Multiple photon packets received by individual 
neighbours are merged in Nc (tessellating) reception cones separately for each frequency u, which limits the number of photon packets 
stored at each particle to at most A^c X A^i^. 



packet merging procedure that respects the chosen angular 
resolution and renders the RT computation time indepen- 
dent of the source number. The merging is done by bin- 
ning photon packets in angle (according to their propagation 
directions) using A^c (tessellating) reception cones. Binned 
photon packets define a single new photon packet per re- 
ception cone whose propagation direction is given by the 
weighted sum^ of the propagation directions of photon pack- 
ets received in that cone. The merging is done separately 



^ We note that the expression in Paper I (Sec. 4.2.3) for the prop- 
agation direction n^^j^ of this new photon packet contains a typo. 
The expression used in that publication, n^^j^ = ^WpUp/ ^Wp, 
where are the unit vectors that represent the propagation di- 
rections of photon packets that are to be merged and Wp are as- 
sociated weights, does generally not result in a unit vector, which 
is inconsistent with the employed notation. However, a unit vec- 
tor representing the propagation direction of the merged photon 
packet can be obtained by an additional explicit normalisation, 

i.e. rim, A; ^m,k/\^m,k\- 



for photon packets of different frequencies. Thus, it does 
not change the mean free path of the photon packets, which 
is important for simulations that contain radiation sources 
with a broad range of spectral properties (e.g., quasars and 
stellar sources). Thanks to the merging, at each particle at 
most Nc X Ni, photons packets need to be stored^. 

^ In principle it is possible to choose the solid angle of the trans- 
mission cones independently of the angular resolution Att/Nc im- 
plied by the emission/reception cone tessellation. The solid angle 
of the transmission cones is the main parameter that determines 
the angular resolution of the photon transport, while the num- 
ber A^c of emission/reception cones controls the number of photon 
packets that need to be stored at each particle. Choosing the solid 
angle of transmission cones smaller than An/Nc would thus result 
in a higher angular resolution while keeping the memory needed 
to store photon packets unchanged. We have successfully tested 
this option by repeating Test 4 in Paper I with Nc = S and a 
transmission cone solid angle of 47r/128. We found the results of 
this simulation to be indistinguishable from the simulation that 
employed A^c = 128 and transmission cone solid angles of 47r/128. 
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The photon transport is performed using RT time steps 
Atr (see Sec. 5.2.3 in Paper I for a detailed discussion of the 
size of time steps in reionisation simulations). During each 
such time step, photons are propagated and their interac- 
tions with the gas are computed until a certain stopping 
criterion is satisfied. The criterion depends on whether one 
aims to solve the time-independent or the time-dependent 
RT equation. In the first case, photons are propagated until 
they are absorbed or have left the computational domain. In 
the second case, photon clocks associated with each photon 
packet are used to synchronise the packet's travel time with 
the simulation time such that photon packets travel at the 
speed of light. 

After each time step, the state of the SPH particles is 
updated according to the interactions (absorptions, scatter- 
ings) with photon packets they experienced. In Sec. 4 we will 
explain, for the example of absorption of ionising radiation, 
how to combine the photon transport discussed here with 
the evaluation of the interactions according to the optical 
depth encountered by photon packets. Finally, the RT time 
is advanced, which concludes the algorithm. A schematic 
summary of the RT method is depicted in the fiow chart in 
Fig. 1. 



3 IONISING PHOTONS IN PRIMORDIAL 
GAS: THEORY 

In this section we outline the physical processes that de- 
termine the evolution of the ionisation state (Sec. 3.1) and 
temperature (Sec. 3.2) of primordial gas exposed to ionising 
radiation. We discuss the underlying equations and present 
the references to atomic data required to evaluate them. The 
description of the numerical implementation used to solve 
these equations is deferred to Sec. 4. 

Readers familiar with the physics of ionisation, recom- 
bination, heating and cooling may wish to skip Sees. 3.1 and 
3.2 and refer to them only when needed. For those readers 
we have summarised the physical processes that we include 
in the computations of the ionisation and thermal state of 
gas in the RT simulations presented later in this work, to- 
gether with the references to the (fits to) atomic data sets 
employed for their numerical evaluation, in Table 1. A de- 
tailed discussion of our choice for certain (fits to) atomic 
data sets and a comparison with other works can be found 
in Pawlik (2009). 

We start with some definitions that we will employ 
throughout the rest of this work. We consider an atomic 
gas of total number density n = ne + ^na^ where ria is the 
number density of species a and Ue is the number density of 
free electrons. The number density ria is related to the total 
mass density p through ria = X^p/ (paTnu), where Xa. is the 
mass fraction of species a and ficx = rria/mu is its mass nia 
in units of the hydrogen mass mu- We assume that the gas 
is of primordial composition, i.e. a G {HI, HII, Hel, Hell, 
Helll} and Xa -h X^e = 1- We wiU set Xa = 0.25. We wih 
make frequent use of the species number density fractions 
with respect to the hydrogen number density, r]a = na/nn 
and the electron fraction r]e = ne/nu- 

In this work, however, we will not make use of this cone decou- 
pling option. 



3.1 Ionisation and recombination 

The evolution of the ionisation state of primordial gas in the 
presence of a photoionising radiation background of mean 
intensity Ju(iy) is determined by the set of rate equations 

drjm 



dt 

dt 
drjneiu 
dt 



= o^HunerjHn - '^Hi(r^Hi + FeHirie) (1) 

= aHeimeryHell — 77Hel(r^HeI + FeHeme) (2) 
= 77HeIl(r^HeII + FeHeime) — Q!HeIime77HeIII, (3) 



supplemented with the closure relations 
Tjm + rjmi = 1 

?7HeI + ?7HeII + ?7HeIII = V^e 
r]mi + ^Hell + 2?7HeIII = r]e, 



(4) 
(5) 
(6) 



where T^a is the photoionisation rate implied by the mean 
intensity of the ionising background and Fe^ and aa are 
the collisional ionisation and recombination rate coefficients 
for species a; rjHe = riHe/riH = ^He(mH/mHe)/(l--^He) de- 
notes the helium abundance (by number); mu and mue are 
the masses of the hydrogen and helium atoms, respectively. 

The ionisation and recombination rates are discussed in 
more detail below. 

3.1.1 Ionisation 

The photoionisation rate T^a determines the number of pho- 
toionisations of species a per unit time and unit volume 
It is defined by (e.g., Osterbrock 1989) 



f 



aiy — ; cr^ 



(7) 



where a G {HI, Hel, Hell}, hp is Planck's constant, cr^a(z^) 
is the photoionisation cross-section and hpUa is the ionisa- 
tion potential. We use the fits to the photoionisation cross- 
section from Verner et al. (1996). Note that hpUui = 13.6 eV, 
hpUKei = 24.6 eV and /ipZ^Heii = 54.4 eV. 

The photoionisation rates can be written as 

4:7rJjy(iy) 



df- 



(8) 



where (a^a) is the average, or grey, photoionisation cross- 
section, 



dv G^a. 



■J Vrv 



dv 



hr, 



.(9) 



We will employ the grey photoionisation cross-section of hy- 
drogen in some of our RT simulations in Sec. 5. For reference 
we note that its value for a blackbody spectrum of temper- 
ature Tbb = 10^ K is (fj^Hi) = 1.63 X 10"^^ cm^ 

In addition to photoionisations we include collisional 
ionisation of HI, Hel and Hell by electron impact. To com- 
pute the corresponding ionisation rates, we employ the fits 
to the collisional ionisation coefficients provided by Theuns 
et al. (1998). 



3.1.2 Recombination 

We write the number of radiative recombinations per unit 
time and unit volume of species a (with a G {HII, Hell, 



© RAS, MNRAS 000, 1-23 



6 A. H. Pawlik and J. Schaye 



Helll}) to energy level I of the recombined species as 

Two radiative recombination coefficients are of special 
interest and are referred to as case A and case B. The case 
A recombination coefficient aAa = ^i=i (^ai is the sum of 
all the recombination coefficients aai. On the other hand, 
the case B recombination coefficient is defined as aBa = 
X^z=2 thus does not include the contribution from 

recombinations to the ground state. 

The introduction of the case B recombination coeffi- 
cient is motivated by the observation that for pure hydrogen 
gas that is optically thick to ionising radiation, recombina- 
tions to the ground state are cancelled by the immediate re- 
absorption of the recombination photon by a neutral atom 
in the vicinity of the recombining ion. RT simulations of ion- 
ising radiation in an optically thick hydrogen-only gas may 
therefore work around the (generally computationally ex- 
pensive) explicit transfer of recombination photons by sim- 
ply employing the case B (instead of the full, i.e. case A) re- 
combination coefficient. Note that this on-the-spot approx- 
imation (e.g, Osterbrock 1989) is only strictly valid when 
considering the transport of ionising radiation in optically 
thick gas, whereas the gas in RT simulations typically shows 
a range of optical depths, from optically thick to optically 
thin. 

To keep the description of our method and its test sim- 
ple and to allow for a detailed comparison with published 
reference results, we will also assume the on-the-spot ap- 
proximation and use case B recombination rates. The ex- 
plicit transport of recombination radiation will be the sub- 
ject of future work. We use the following coefficients to 
describe radiative recombinations (Table 1). For HII and 
Helll radiative recombination, we employ the fits from Hui 
& Gnedin (1997). For the Hell radiative recombination co- 
efficient, we employ the tabulated coefficients of Hummer & 
Storey (1998) using linear interpolation in log- log. 

We have not yet discussed the dielectronic contribution 
to the HeH recombination coefficient. Dielectronic HeH re- 
combination (e.g. Savin 2000; Badnell 2001 for a review) 
is the dominant recombination process for temperatures 
T > 10^ K. We therefore add the dielectronic contribution 
to the HeH recombination rates, making use of the fit pre- 
sented in Aldrovandi & Pequignot (1973). 



3.2 Heating and cooling 

Our main goal in this work is to present and test an imple- 
mentation of TRAPHIC to compute, in addition to the ionisa- 
tion state, the evolution of the temperature of gas exposed 
to ionising radiation. For the discussion it is helpful to re- 
view the relevant thermodynamical relations, which is the 
subject of this section. 

The internal energy per unit mass for gas composed of 
monoatomic species that are at the same temperature T is 



2 p 



3 /cbT 
2 firnn ' 



(10) 



where /cb is the Boltzmann constant and /i = p/{nmu) is 
the mean particle mass in units of the hydrogen mass. 
From the first law of thermodynamics 



where P is the pressure, V the volume and 1-L and C are 
the radiative heating and radiative cooling rates, normalised 
such that the rates of energy gain and loss per unit volume 
are described by n^Ti and n^C, respectively. Assuming that 
d(pV) = 0, as is the case for an SPH particle, it follows that 



dt~ pV dt ^ p 



3.2.1 Cooling 



(12) 



The normalised cooling rate C is the sum of the normalised 
rates of the individual radiative cooling processes. We in- 
clude all standard cooling processes: collisional ionisation by 
electron impact, radiative and dielectronic recombination, 
collisional excitation by electron impact, bremsstrahlung 
and Compton scattering. The expressions for the corre- 
sponding cooling rates are taken from the references listed 
in Table 1. 



3.2.2 Heating 

The normalised heating rate Ti is the sum of the rates of the 
individual radiative heating processes. In the following we 
only consider the contribution from photoionisation heating, 
which will be the main contributor to the heating rate for 
the high-redshift RT simulations of interest. We, however, 
note that Compton heating by X-rays may not be negligible 
(Madau k Efstathiou 1999). 

We write the heating rate due to photoionisation as 



^H^7 — (^HI^7HI + ^HeI^7HeI + ^HeII^7HeIl)^H , 

where 

47rJz.(zy) 



8^ 



- r 



dv- 



(13) 



(14) 



is the heating rate per particle of species a. Using Eq. 7, we 
can write 



^7a — F-yo; ( 

where 



(15) 



[f 

.J Vol 

[/° 



, 47rJ^(z/) 

dv ; —^Gry 



t{y)(hpV — hpi^a) 



dv — (j^<x{y) 



(16) 



is the average excess energy of absorbed ionising photons. 
For reference, the average excess energy for photoionisation 
of hydrogen, assuming a blackbody spectrum of temperature 
Tbb = 10^ K, is (em) = 6.32 eV. 

Sometimes, e.g. when considering the energy balance of 
entire HH regions, one is interested in the total photoheat- 
ing rate integrated over a finite volume, assuming all photons 
entering this volume are absorbed within it. The average ex- 
cess energy injected at each photoionisation in this optically 
thick limit is also obtained from Eq. 16, but after setting 
C7q;(z^) = 1, since all photons are absorbed (e.g., Spitzer 
1978, p.135), 

/ thick\ 



d{upV) = -PdV -f- n|CH - C)V, 



(11) 



[/: 
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dv 



dv 
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hpV 



(17) 
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For reference, the value of the average excess energy for pho- 
toionisation of hydrogen in the opticaUy thick hmit, assum- 
ing a blackbody spectrum of temperature Tbb = 10^ K, is 

In writing Eqs. 16 and 17 we assumed that ah of the 
photon excess energy is used to heat the gas, corresponding 
to a complete thermalization of the electron kinetic energy. 
In reality, (very energetic) photo-electrons may lose some of 
their energy due to the generation of secondary electrons 
(e.g., Shull k, van Steenberg 1985; Furlanetto k, Stoever 
2010). 



4 IONISING PHOTONS IN PRIMORDIAL 
GAS: IMPLEMENTATION 

Here we extend the description of traphic given in Sec. 2 
to the transport of ionising photons by describing our imple- 
mentation of the absorption of ionising photons and the sub- 
sequent computation of the species fractions and gas tem- 
peratures. 



4.1 Absorption of ionising photons 

The number of ionising photons that are absorbed during 
the propagation of a photon packet over distance dij be- 
tween neighbouring particles i and j is given by SNg,hs,iy = 
SNin,u[l — exp(— t(z/))], where 6Nin,u is the number of pho- 
tons inside the photon packet before propagation and the 
optical depth r(iy) is the sum t(i/) = ^a^cx(i^) of the op- 
tical depths of each absorbing species a G {HI, Hel, Hell} 
and 



i: 



dr acx(iy)ncx{r) 



^{-u)nc,{Yj)di. 



(18) 



The last approximation is reasonable because SPH neigh- 
bours will have similar densities. We consider these photons 
to be absorbed by particle j. 

The number of photons (5A^abs,a(i^) absorbed by species 
a is determined from the number of absorbed photons 
SN 3,hs,u using 



(19) 



where we choose (un-normalized) weights Wa(iy) = Ta{iy) 
(Osterbrock 1989; see also, e.g., Trac & Cen 2007). The to- 
tal number of photons AA^abs,c^(z^) absorbed by species a is 
the sum of the photons absorbed due to the propagation of 
photon packets from all neighbouring particles during the 
RT time step Atr, i.e., AA^abs,a(i^) = J2^^a.hs,cx{T^)- 

We pause to note that other choices for the weights 
Wa (as employed in, e.g., Bolton, Meiksin, & White 2004; 
Maselh, Ferrara, & Ciardi 2003; Whalen & Norman 2008) 
will in general imply an unphysical distribution of the num- 
ber SNohs.u of absorbed photons amongst the individual 
species. To see this, consider, for example, a gas parcel with 
optical depth r that contains only hydrogen atoms. Let an 
arbitrary fraction /a of the hydrogen atoms be labelled A 
and let the remaining fraction = 1 — of the hydrogen 
atoms be labelled B. Suppose that the optical depth of the 
hydrogen atoms with label A is ta- The ratio of the number 
of photons absorbed by hydrogen atoms with label A to the 



number of photons absorbed by all hydrogen atoms should 

be 



(5iVabs,a(z^) _ Ta(z^) 



r(^) 



while Eq. 19 implies 



^iVabs,: 



(20) 



(21) 



The two ratios generally agree only if oc , where a takes 
values A and B. Other choices of the weights Wa would im- 
ply that the probability for absorption of ionising photons 
by hydrogen atoms with label A is different from the proba- 
bility for absorption of ionising photons by hydrogen atoms 
with label B. This is unphysical, since there is no physical 
difference between these two types of hydrogen atoms. 

As described in Paper I, photons absorbed by a virtual 
particle (ViP) are redistributed amongst the iVngb neigh- 
bouring SPH particles that have been used to compute its 
species densities. This is necessary, because ViPs are tem- 
porary constructs; physical properties like the gas species 
fractions are only defined and stored for the SPH particles. 
There, however, is an important change with respect to the 
original description. Previously, we assigned to each of the 
neighbours a fraction of the absorbed photons that is pro- 
portional to the value of the ViP's SPH kernel W at its 
position. In the current version we assign, to each of the 
neighbours, a fraction of the photons absorbed by species a 
that is proportional to the neighbour's contribution to the 
SPH estimate of the density of that species at the location 
of the considered ViP. The current version is equivalent to 
the original version of Paper 1 ii X = 1 (i.e. no helium) 
and if all neighbours have the same neutral hydrogen mass. 
However, in general this will not be true, in which case the 
current version is the only self- consistent one. We discuss 
the differences between the current and the original version 
in detail in App. A. 

The number of photons ANahs,a{i^) that are absorbed 
by a given particle during the RT time step Atr is used to 
obtain the photoionisation rates Fa for that particle directly 
(i.e., without reference to the mean intensity Ji^) using 



rjcyAfnAtrTjcy = AA/'abs,a(^), 



(22) 



where A/h = mXu/mu is the number of hydrogen atoms 
associated with the SPH particle of mass m. The photoion- 
isation rates are then used to advance the species fractions 
and the gas temperature as we will describe in the next sec- 
tion. 



4.2 Integration of the rate equations 

Here we present our numerical method to solve the equations 
of the evolution of the ionisation balance and temperature 
of gas exposed to ionising radiation (Eqs. 1-6 and Eq. 12). 
This method is an extension of the subcycling method de- 
scribed in Paper I, which we therefore briefly recall. 

In Paper I we presented a method to follow the ion- 
isation state of a (hydrogen-only) gas parcel exposed to 
(hydrogen-) ionising radiation at fixed temperature. The 
ionisation rate equations were solved at the end of each 
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Table 1. References to (fits to) the atomic data used to calculate photoionisation rates, collisional ionisation rates, recombination rates 
and cooling rates in the simulations presented in this work. See Pawlik (2009) for detailed discussions and comparisons of our choices in 
favour of certain references. 



Photoionisation 
Collisional ionisation 


HI, Hel, Hell photoionisation cross-sections 

HI, Hel, Hell collisional ionisation rate coefficients 


Verner et al. (1996) 
Theuns et al. (1998) 


Recombination 


HII, Helll recombination rate coefficients 

Hell recombination rate coefficient 

Hell dielectronic recombination rate coefficient 


Hui & Gnedin (1997) 
Hummer & Storey (1998) 
Aldrovandi & Pequignot (1973) 



Collisional ionisation cooling 
Collisional excitation cooling 
Recombination cooling 



Cooling by bremsstrahlung 
Compton cooling 



HI, Hel, Helll collisional ionisation cooling rate 
HI, Hel, Helll collisional excitation cooling rate 
HII, Helll recombination cooling rate (case A and B) 
Hell recombination cooling rate (case A and B) 
Hell dielectronic recombination cooling rate 
Bremsstrahlung cooling rate 
Compton cooling rate 



Shapiro & Kang (1987) 

Cen (1992) 

Hui & Gnedin (1997) 

Hummer & Storey (1998) 

Black (1981) 

Theuns et al. (1998) 

Theuns et al. (1998) 



RT time step Atr by explicit numerical integration (here- 
after also referred to as subcycling) using subcycle steps 
St = min(/req, Atr), where 



is the time scale to reach ionisation equilibrium (Eq. 20 in 
Paper I), Tree = l/(^e<^Hii) IS the recombination time scale, 
Tion = l/(r^Hi + TieTeHi) IS the ionisatiou time scale and / 
is a dimensionless factor that controls the integration accu- 
racy. Subcycling allows the RT time step Atr to be chosen 
independently of the values of the ionisation and recombi- 
nation time scales on which the species fractions evolve. A 
RT time step Atr limited by the ionisation and recombina- 
tion time scales would prevent efficient RT simulations since 
these time scales may become very small. 

In this work we are interested in the self-consistent com- 
putation of the non-equilibrium ionisation state of gas with 
an evolving temperature. As for the case of a non-evolving 
temperature studied in Paper I, we integrate the ionisa- 
tion rate equations over subcycle steps 6t = min(/req, Atr). 
The time scale Teq to reach ionisation equilibrium is com- 
puted using Eq. 23 with Tion = 0,(^70^ + '^e^ea) and 
Tree = 1/ ^eO^a- Rccombiuation and collisional ionisation 
rates are determined using the temperature at the beginning 
of each subcycle step and the species fractions are advanced 
in a photon-conserving manner as detailed^ in Paper I. 

In addition, the temperature is advanced by evolving 
the internal energy according to Eq. 12 over the same sub- 
cycle step assuming isochoric evolution {dV = 0), which is 
appropriate for a fixed gas distribution (and thus during a 
single hydro-step in radiation- hydrodynamical simulations). 
We use the mean particle mass /i derived from the current 
species fractions to convert between temperature (which is 
required to compute the rate coefficients) and internal en- 
ergy using Eq. 10. Note that the species fractions and the 
temperature are evolved independently of each other over 
a single subcycle step of size 6t. We thus implicitly assume 
that during any of these steps the species fractions and the 

^ In Paper I we only considered ionisation of gas of pure hydro- 
gen. The corresponding expressions, however, are straightforward 
to generalise to include the ionisation of helium. 



temperature do not evolve significantly. This assumption is 
excellent because the species fractions and the temperature 
evolve, by definition, on time scales large compared with the 
size of the subcycle steps. The evolutions of the species frac- 
tions and the temperature are coupled at the beginning of 
the next subcycle step, where the new species fractions and 
the new temperature determine new collisional ionisation, 
recombination and cooling rates. 

We now describe our numerical implementation of the 
subcycling. We limit ourselves to the description of how we 
advance the internal energy over a single subcycle step as 
the implementation of the subcycling of the species frac- 
tions was already described^ in Paper I. The internal energy 
is advanced by solving a discretized version of the energy 
equation (i.e., Eq. 12 with dV = 0). We make use of im- 
plicit Euler integration when the subcycle step is larger than 
the time scale Tu = u/{du/dt) on which the internal energy 
evolves. That is, if St > Tu we advance the internal energy 
according to 

ut+5t = ut + -^{l-Lt+5t - Ct+5t)St. (24) 
Pt 

The last equation is solved iteratively for Ut+5t by finding 
the zero of the function 

f(ut+5t) = Ut+5t -Ut -{l-it+5t - Ct+5t)St (25) 

Pt 

and setting l-Lt+5t — 1-it and assuming Ct-\-5t = Ct during the 
first iteration. If, instead, St < Tu, we employ the explicit 
Euler integration scheme, 

2 

ut+st = + — (Ht - Ct)St, (26) 
pt 

Our implementation combines the advantages of the explicit 
scheme (its accuracy) with that of the implicit scheme (its 
stability; see, e.g., Shampine & Gear 1979 and Press et al. 
1992 for useful discussions on implicit and explicit integra- 
tion). 

In Paper I (for the case of a constant temperature), we 
sped up the subcycling of the species fractions by keeping 
the species fractions fixed once ionisation equilibrium has 

^ See footnote 8 
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been reached . We employ a similar recipe here. However, 
thermal equilibrium is reached on the time scale Tu, which 
may be much larger than the time scale Teq to reach ionisa- 
tion equilibrium. In this case the temperature continues to 
evolve after the species fractions have attained their (quasi-) 
equilibrium values. The evolution of the temperature implies 
an evolution of the recombination and collisional ionisation 
rates, and hence an evolution of the equilibrium ionisation 
balance. Our recipe for speeding up the subcycling should 
respect this evolution. 

We therefore proceed as follows. Once ionisation equilib- 
rium has been reached, we stop the subcycling of the species 
fractions. Over the remainder of the time step Atr only the 
internal energy is subcycled, which can be done using time 
steps = /u X Tu, where /u < 1 is a dimensionless param- 
eter that controls the accuracy of the integration (we set 
fu — f)' This results in a speed-up since typically Sut ^ St. 
After each such subcycle step, we reset the species fractions 
to their current equilibrium values. The equilibrium species 
fractions are obtained by iteratively solving the set of equa- 
tions 1-6 with drja/dt = 0. 

In summary, we solve the evolution of the ionisation 
balance and temperature using a hybrid numerical method 
that makes use of both explicit and implicit Euler integra- 
tion schemes. The ionisation rate equation is solved explic- 
itly using the subcycling procedure presented in Paper I. 
This ensures the accurate conservation of photons and al- 
lows one to choose the size of the RT time step indepen- 
dently of the (often very small) ionisation and recombination 
time scales, a pre-requisite for efficient RT simulations. The 
temperature is evolved along with the ionisation balance by 
following the evolution of the internal energy of the gas. We 
use an explicit discretisation scheme to advance the inter- 
nal energy if the cooling time is larger than the size of the 
subcycle step. For smaller cooling times, stability consider- 
ations lead us to employ an implicit discretisation scheme 
to advance the internal energy. Once ionisation equilibrium 
has been reached, the subcycling computation is sped up by 
fixing the species fractions to their (temperature-dependent) 
quasi-equilibrium values. From then on, only the evolution 
of the internal energy is subcycled. 



5 IONISING PHOTONS IN PRIMORDIAL 
GAS: RESULTS 

In this section we perform simulations to test our new, ther- 
mally coupled implementation of TRAPHIC. We also use these 
simulations to discuss the differences in RT simulations per- 
formed using the grey approximation in the optical thick 
and thin limits and compare simulations using the grey ap- 
proximation to simulations that solve the RT using multiple 
frequency bins. 

We begin in Sec. 5.1 with verifying that the subcy- 
cling method described in Sec. 4.2 can be successfully em- 
ployed to solve for the non-equilibrium ionisation balance 
and temperature of gas exposed to ionising radiation. Then, 
in Sec. 5.2, we present a set of reference solutions for the 

We assume that ionisation equilibrium is reached once either 
the fractional change in all species fractions individually becomes 
smaller than a predefined small value (here we use 10~^). 



idealised problem of a single spherically symmetric expand- 
ing HII region that we will employ later to test the perfor- 
mance of TRAPHIC in this same problem. We compare ref- 
erence solutions derived in the grey approximation with the 
exact mult i- frequency results and discuss their differences. 
Thereafter, in Sees. 5.3 and 5.4, we investigate traphic's 
performance in RT simulations of increasing complexity: in 
Sec. 5.3 we compute the ionised fractions and temperatures 
around a single source in a homogeneous density field and in 
Sec. 5.4 we follow the ionising radiation of multiple sources 
in a highly inhomogeneous density field. Throughout we will 
compare the results obtained with TRAPHIC to analytical and 
numerical reference results. 

The simulations were performed with TRAPHIC imple- 
mented in a modified version of GADGET-2 (Springel 2005). 
All simulations were run on static density fields. We also 
remind the reader that, to facilitate comparisons with refer- 
ence simulations, we do not explicitly follow recombination 
radiation but treat it using the on-the-spot-approximation. 



5.1 Test 1: Subcycling 

Here we test the subcycling approach to the computation 
of the coupled evolution of the non-equilibrium ionisation 
balance and temperature of gas exposed to ionising radi- 
ation that we have introduced in Sec. 4.2. Our aim is to 
demonstrate that, given a fiux impinging on a gas parcel 
(or, equivalent ly, a photoionisation rate experienced by this 
parcel), the subcycling allows for an accurate computation 
of the evolution of its ionisation state and temperature, in- 
dependent of the size of the RT time step Atr- 

The setup of the test is as follows. We simulate the evo- 
lution of the ionisation state of an optically thin gas parcel 
with hydrogen number density nn = 1 cm~^. For simplic- 
ity and clarity of the presentation, we set the hydrogen mass 
fraction to X = 1 (i.e., no helium). The simulation starts at 
time t — with a fully neutral parcel with initial temper- 
ature T — 10^ K. We then apply a photoionising flux of 
F = 10^^ s~^ cm~^ with a blackbody spectrum of char- 
acteristic temperature Tbb = 10^ K. Consequently, the par- 
cel becomes highly ionised and is heated to a temperature 
T - 10^ K. After t = 50 Myr we switch off the ionising flux 
and the parcel recombines and cools. The simulation ends at 
^end = 1 Gyr. The test here is identical to Test presented 
in Iliev et al. (2006a), except for the switch-off time (Iliev 
et al. 2006a used tend = 0.5 Myr). 

We employ a grey photoionisation cross-section 
(cr^Hi) = 1-63 X 10~^^ cm^ (Sec. 3.1), yielding a photoion- 
isation rate F^hi = 1.63 x 10~^ s~^. We assume that each 
photoionisation adds (chi) = 6.32 eV to the internal energy 
of the gas (Sec. 3.2), which corresponds to the optically thin 
limit. We solve the equations for the evolution of the ionisa- 
tion state and temperature of the gas parcel by subcycling 
them on subcycle time steps 5t over consecutive time inter- 
vals At. Note that in a full RT computation these intervals 
would correspond to the RT time steps Atr. Here we distin- 
guish between At and Atr only because in this test we are 
considering a single gas parcel with prescribed photoionisa- 
tion rate and do not perform RT simulations. The dimen- 
sionless parameter / that controls the size of the subcycling 
steps 8t (and hence the accuracy of the subcycling) is set to 
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Figure 2. Test 1. Single gas parcel with prescribed photoionisation rate. The top panels show the evolution of the neutral hydrogen 
fraction (left) and temperature (right) for simulations with different time steps At (in the range 5 x 10~^ yr — 5 x 10^ yr), as indicated 
in the legend in the top left panel. All simulations start at time t = and evolve the gas parcel on subcycle steps 6t = min(/req, At), 
where Teq is the time scale on which the gas parcel approaches equilibrium and / = 0.01 is a dimensionless parameter that controls the 
accuracy of the subcycling. Note that the first output occurs only at the end of the first time step, i.e., at time t = At, which is different 
for the different simulations. The bottom panels show the relative errors of the evolution shown in the top panels with respect to the 
evolution obtained from the simulations with the next smaller time step. The relative errors are small, < 1%, except during a brief phase 
of rapid evolution in the ionised fraction. The relative error can be reduced by choosing a smaller value for the accuracy parameter /. 



/ = 10 ^. When computing Compton cooling rates off the 
cosmic microwave background, we assume a redshift z = 0. 

Fig. 2 shows the results. The top left and top right 
panels show the evolution of the neutral hydrogen fraction 
and the temperature, respectively, for simulations with time 



steps At = 5 X (10" 



10- 



10-\ 10\ 10^ 10^) yr. Note 



that in order to limit the computation time, not all of the 
simulations have been evolved until the end of the simu- 
lation time. They were stopped once their simulation time 
overlapped with that of the simulation with the next larger 
time step. The earliest output of a given simulation is at 
t = At (and hence, in Fig. 2, different curves start at dif- 
ferent times), but all simulations started at t = and were 
numerically evolved on subcycle time steps St < At as de- 
scribed in Sec. 4.2. 

The gas parcel quickly approaches photoionisation equi- 
librium, reaching its equilibrium neutral fraction after ^10 



(photo-) ionisat ion time scales (ric 



= r 



7HI 



0.02 yr). 



During this period, photoheating raises its temperature to 
T 2 X 10^ K. Around t ~ 10^ yr, the neutral fraction ex- 
hibits a slight decrease. As noted in Iliev et al. (2006a), this 
behaviour is caused by the decrease in the recombination 
rate due to the rise in temperature that can be observed at 
this time. The fact that the temperature still evolves after 
the neutral fraction reached its equilibrium value means that 
thermal equilibrium is reached on a larger time scale than 
photoionisation equilibrium. Note that the thermal equilib- 
rium phase is missed in the test simulations of Iliev et al. 
(2006a) because these simulations were stopped at a much 
earlier time. 

The observed behaviour can be understood as follows. 
When thermal equilibrium is approached from a tempera- 
ture lower than the equilibrium temperature, the net cool- 
ing rate is approximately given by the photoheating rate. 
In photoionisation equilibrium, the photoheating rate is 
proportional to the recombination rate. The time scale 



Tu = u/(du/dt) to reach thermal equilibrium can there- 
fore be expressed in terms of the recombination time Tree = 

l/(neaHIl), 



Tu 



(3/2)n/cBT 

(3/2)nkBT 
nmineami{em) 
(3/2)n/cBT 

/ \ Tree 

{em/nmi 

Tree 5 



(27) 

(28) 

(29) 
(30) 



where in the last step we assumed that the gas is highly 
ionised, i.e. riHii ^ nu ^ n/2, and that T ^ 10^ K. 
The recombination time (and hence the cooling time) is 
much larger than the time Teq = (TionTree)/(rion + Tree) to 
reach ionisation equilibrium which asymptotes to Tion for 
Tion ^ Tree (scc also the discussiou in Sec. 5.1 in Paper I). 
Here, Tion = ^ ^-^^ "^^ee ~ 10^ yr. Accordingly, 

thermal equilibrium is reached much later than photoioni- 
sation equilibrium. 

After thermal equilibrium has been reached, the ionis- 
ing flux is switched off and the parcel recombines and cools. 
Once it has cooled to a temperature T < 10^ K, cooling be- 
comes inefficient. The temperature of the recombining parcel 
therefore remains roughly constant. 

In the bottom panels of Fig. 2 we quantify the accuracy 
of our subcycling approach. Ideally, we would like to com- 
pare the numerical results to an exact analytical reference 
solution. However, such a solution exists only for the case 
of a constant temperature^^ (see, e.g., the appendix in Dove 

We mention that by repeating the test at fixed temperature, 
we have convinced ourselves that the ionisation history computed 
using our subcycling recipe follows the analytical solution very 
closely (see Pawlik 2009). 
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& Shull 1994). Instead, we therefore show the relative er- 
ror of the evolutions shown in the top panel with respect to 
the evolutions obtained from the simulation with the next 
smaller time step. 

For all our choices of the time step At and for most of 
the simulation time the relative errors are small, < 1% or 
even much smaller. During the initial phase of rapid evolu- 
tion the relative error in the ionised fraction briefly becomes 
as large as 10%. In practice, such errors will have little im- 
pact on the results of RT simulations if photons are con- 
served and if the equilibrium solution is still obtained with 
high accuracy (as is the case). Moreover, relative differences 
of order 10% are already implied by uncertainties in current 
atomic data used to compute the ionisation and recombina- 
tion rates and the radiative heating and cooling rates (as will 
be discussed in Sec. 5.2.2). Note that the relative error can 
be reduced by lowering the numerical factor /, which con- 
trols the size of the subcycle steps and hence the integration 
accuracy. We conclude that the results of the subcycling are 
insensitive to the size of the simulation time step. 

In summary, we have demonstrated that our subcycling 
recipe accurately computes, independently of the size of the 
RT time step, the combined evolution of the neutral frac- 
tion and temperature of gas exposed to hydrogen-ionising 
radiation. In the following sections we will employ the sub- 
cycling to compute the species fractions and temperature of 
gas particles in RT simulations. 



5.2 HII region expansion. Reference results and 
comparisons of multi-frequency and grey 
solutions 

In the next section (Sec. 5.3) we will apply our new im- 
plementation of TRAPHIC to compute the evolution of the 
ionisation state and temperature around an ionising source 
surrounded by gas of constant density. This is an idealised 
test problem designed to facilitate the verification of our 
implementation through the direct comparison to results 
obtained with an improved version of our one-dimensional 
(1-d) RT code (Pawlik k Schaye 2008; hereafter referred 
to as ttId, which stands for TestTraphiclD), which solves 
the rate equations using the same techniques (and code) as 
TRAPHIC, as well as to published reference results obtained 
with other RT codes for the same test problem (Iliev et al. 
2006a). In this section we present these reference results. We 
also discuss the applicability of the grey approximation for 
solving mult i- frequency RT problems. 

We start in Sec. 5.2.1 by presenting reference solu- 
tions obtained with our 1-d RT code ttId for the case 
of hydrogen-only gas at fixed temperature. This is an im- 
portant case because it allows analytical solutions to be 
derived against which the numerical results obtained with 
ttId can be compared. Then, in Sec. 5.2.2, we compare the 
performance of ttId in a simulation of hydrogen-only gas 
in which the gas temperature is allowed to evolve due to 
photoheating and radiative cooling to published numerical 
results obtained for the same problem. Finally, in Sec. 5.2.3, 
we discuss results from simulations with ttId in which the 
gas also contains helium. 



5. 2. 1 HII region in pure hydrogen gas at fixed temperature 

In this section we discuss the RT problem of an expand- 
ing HII region. Despite its simplicity, an analytical solution 
to this problem cannot generally be obtained, even if the 
gas densities are assumed to be non-evolving (as is the case 
throughout this work). This is because the coupling between 
the ionisation and temperature state through the depen- 
dence of the collisional ionisation, recombination and cool- 
ing rates on the temperature and species fractions impedes 
the evaluation of the governing differential equations (Eqs. 1 
- 6 and 12). 

To provide an approximate point of reference, we re- 
call the evolution of an HII region at fixed gas temperature, 
for which an analytical solution is known (under the ap- 
proximation that the ionised region is fully ionised; we will 
also ignore collisional ionisations, although this is not neces- 
sary). We have reviewed this solution in Paper I, where we 
showed that the radius of the ionised sphere around a source 
of ionising luminosity that is located in a homogeneous 
hydrogen-only medium of density nn is given by 

rI(^) = rs(l-e-*/"=)l/^ (31) 

where rg = [3A/'7/(Q!BHimH)]^'^^ is the Stromgren radius and 
Ts = l/(aBHimH) is the Stromgren time scale, which equals 
the recombination time for fully ionised gas. In some of our 
comparisons we will employ this approximate point of ref- 
erence. We will refer to it as an analytical approximation. 
On the other hand, because of the lack of an accurate ana- 
lytical solution, we will mostly employ results obtained with 
our 1-d RT code ttId in our benchmarking below. For this 
reason, we will first discuss its performance. 

We start by verifying our multi-frequency treatment in 
ttId by comparing its performance in a simple HII region 
test problem to the corresponding equilibrium solution that 
can be derived analytically (except for a numerical evalua- 
tion of the integrals involved). The test consists of simulat- 
ing the spherically symmetric growth of the ionised region 
around a single ionising source in a homogeneous hydrogen- 
only medium at fixed temperature. The source has a black- 
body spectrum with temperature 10^ K and emits radiation 
with an ionising luminosity = 5 x 10^^ photons s~^. The 
gas density is nn = 10 ~^ cm~^. The initial ionised frac- 
tion is assumed to be ?7hii — 0, and we use a recombination 
coefficient aBHii = 2.59 x 10~^^ cm^ s~^, independent of 
radius and time. Collisional ionisation is not included. For 
reference, with the physical parameters mentioned above, 
the Stromgren time is rs = 122.4 Myr and the Stromgren 
radius is rs = 5.4 kpc. The spatial resolution, the time step 
and the number of frequency bins used in the simulation 
with ttId are chosen such as to achieve numerical conver- 
gence. 

In Fig. 3 we show the neutral (ionised) fraction profile 
in photoionisation equilibrium. Diamonds show the result of 
the simulation with ttId (at t — 2000 Myr). The blue sohd 
curve indicates the exact equilibrium solution obtained by 
solving (e.g., Osterbrock 1989) 

^HI,eq(r)nH f i xr { \ -^u 2 ( \ 2 /QO^ 

/ diy N^{iy)e "cr^, = ryHii,eq(^)^H<^BHii, (32) 

where the frequency-dependent optical depth Ti^{r) is given 
by 
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Figure 3. Photoionisation equilibrium profiles of the neutral and 
ionised fraction around a single blackbody source in a homoge- 
neous hydrogen-only medium. The numerical result obtained with 
ttId (diamonds) shows excellent agreement with the analytically 
computed exact results (Eq. 32; blue solid curve). For comparison, 
we also show the analytically computed exact solutions assuming 
the grey approximation (red dashed curve) and the monochro- 
matic treatment used in Test 1 of Paper I (black dotted curve). 
The grey approximation agrees with the exact solution in the op- 
tically thin limit (i.e., in the absence of spectral hardening), while 
the monochromatic treatment fails at all distances because it im- 
plies photoionisation rates inconsistent with the spectrum of the 
considered blackbody source. 



Tu{r) = nuCTu / dr r]m,eci{r). (33) 

The simulation result is in excellent agreement with the 
exact equilibrium solution, verifying our mult i- frequency im- 
plementation of ttId. For comparison, we also show the 
exact equilibrium solution assuming that the radiation is 
monochromatic (dotted black curve) with a photoionisa- 
tion cross-section evaluated at the ionisation threshold, i.e. 
C7HI = 6.3 X 10~^^ cm^. We also show the exact equilib- 
rium solutions in the grey treatment, i.e. using the average 
cross-section (ct^hi) = 1.63 x 10~^^ cm^ (dashed red curve). 

The reason for the differences between the results of 
the multi-frequency computation and the results of the 
grey and monochromatic computation can be readily un- 
derstood. The absorption cross-section for ionising photons 
is a strongly decreasing function of the photon energy. The 
ionising photons with the lowest energy are therefore pref- 
erentially absorbed, which leads to an increase in the typi- 
cal photon energy with distance. This effect is referred to as 
spectral hardening. Because the photon mean free path is in- 
versely proportional to the absorption cross-section, spectral 
hardening increases the width of the ionisation front with re- 
spect to that obtained in the absence of spectral hardening. 
Note that spectral hardening only becomes important for 
large optical depths, which explains why the grey approxi- 
mation reproduces the mult i- frequency solution at small dis- 
tances where the optical depth is low. The monochromatic 



approximation, on the other hand, implies an inappropri- 
ate value for the photoionisation rate and hence fails to de- 
scribe the present multi-frequency problem at all distances 
from the source. Note that both the grey and the monochro- 
matic approximation will provide a better description of the 
multi-frequency problem for sources with a softer radiation 
spectrum. 

5.2.2 HII region in pure hydrogen gas with an evolving 
temperature 

Having demonstrated the validity of our multi-frequency 
treatment with ttId, we now repeat the test problem 
from the previous section but this time we account for 
the self-consistent evolution of the gas temperature due 
to photoheating and radiative cooling. The physical pa- 
rameters for the test are taken from Iliev et al. (2006a). 
We consider an ionising source embedded in a homoge- 
neous hydrogen-only density field with number density 
riH — 10~^ cm~^. The source emits jVj = f dv JVj{iy) — 
5 X 10^^ ionising photons s~^ with a blackbody spectral 
shape J\f^{iy) corresponding to a blackbody temperature 
Tbb = 10^ K. The test described here is identical to Test 1 
in Paper I, except that now the gas temperature is allowed 
to vary due to heating and cooUng processes as described in 
Sec. 3.2 (with Compton cooling off the redshift z = cosmic 
microwave background included) and that collisional ionisa- 
tion is included. The gas is assumed to have an initial ionised 
fraction ?7hii = 1.2 x 10~^ (approximately corresponding to 
the ionised fraction implied by collisional ionisation equilib- 
rium at the temperature T = 10^ K). Its initial temperature 
is set to 10^ K. As before, the spatial resolution, the time 
step and the number of frequency bins used in the simu- 
lation with ttId are chosen such as to achieve numerical 
convergence. 

Fig. 4 shows the neutral (ionised) fraction and temper- 
ature profiles at time t = 100 Myr using ttId (black solid 
curves). We compare these multi- frequency results to results 
obtained using the grey approximation, i.e. using an average 
cross-section (ct^hi) = 1.63 x 10""^^ cm^ and grey photoheat- 
ing rates. We employ grey photoheating rates computed in 
the optically thin limit (red dotted curves), according to 
which each photoionisation adds (chi) = 6.32 eV to the in- 
ternal energy of the gas, and in the optically thick limit 
(blue dashed curves), in which case each photoionisation 
adds (cm^^) — 16.01 eV to the internal energy of the gas 
(Sec. 3.2). We employ the labels grey thin and grey thick 
to distinguish the two grey simulations from each other. We 
also show the results obtained (in three-dimensional RT sim- 
ulations) with the RT codes C^-RAY (Mellema et al. 2006), 
CRASH (Ciardi et al. 2001; MaselU, Ferrara, & Ciardi 2003) 
and FTTE (Razoumov & Cardall 2005) for the same test 
problem, as published in Iliev et al. (2006a). 

The differences in the neutral fractions between the grey 
and the multi-frequency simulations that we have discussed 
above for Fig. 3 are again clearly visible (top panel of Fig. 4). 
The grey thin simulation yields results that asymptote to 
those obtained in the multi-frequency simulation at small 
distances from the ionising source. At large distances, i.e. 
near the ionisation front and beyond, on the other hand, 
the multi-frequency simulation implies significantly larger 
ionised fractions than those implied by this grey simulation. 
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Figure 4. Test 2 (using ttId). Comparison of the grey approxi- 
mations with the full multi-frequency solution. The figure shows 
spherically averaged profiles of neutral (ionised) fraction (top) 
and temperature (bottom) at time t = 100 Myr. The black solid 
curve shows the multi-frequency solution. The red dotted (blue 
dashed) curve shows its grey approximation assuming photoheat- 
ing rates computed in the optically thin (thick) limit. The grey 
and multi-frequency simulations show clear differences close to 
and beyond the ionisation front, where the large optical depth 
causes spectral hardening of the emitted blackbody spectrum. 
For reference, we also show results obtained with other RT codes 
as published in Iliev et al. (2006a). The green dash-dotted and 
black solid curves in the top panel are difficult to distinguish 
because they fall nearly on top of each other. The differences 
between these results at large distances mainly reffect the differ- 
ences in the numerical treatment of multi-frequency radiation in 
these codes. Most of the differences close to the ionising source 
have their origin in the use of different assumptions for comput- 
ing photoheating rates, as the comparison to the results obtained 
with ttId reveals. 



This is because the photon mean free path is larger in the 
multi-frequency simulation than in the grey simulations due 
to spectral hardening, leading to a smoother transition of 
the neutral fraction between the highly ionised gas interior 
to and the neutral gas far ahead of the ionisation front. 

The grey thick simulation yields neutral fractions that 
are very similar to those found in the grey thin simulation. 
However, the grey thick simulation yields slightly lower neu- 
tral fractions than the grey thin simulation, since it yields 
slightly larger temperatures, and thus smaller recombina- 
tion rates, throughout the ionised region (bottom panel of 
Fig. 4). In contrast to the grey thin simulation, the neu- 
tral fractions obtained in the grey thick simulation therefore 
do not asymptote to those obtained in the multi-frequency 
simulation at small distances to the ionising source. Instead, 
they remain systematically too small. 

The differences between the grey and multi-frequency 
simulations (and between the grey thin and grey thick 
simulations) become particularly apparent when inspect- 



ing the corresponding temperature profiles. The multi- 
frequency simulation yields substantially higher gas temper- 
atures ahead of the ionisation front. This pre-heating is a 
simple consequence of the increase in the photon mean free 
path above that in the grey simulations. As already noted, 
at fixed radii the grey thick simulation shows systematically 
higher gas temperatures than the grey thin simulation. The 
reason is that in the optically thin limit the contribution 
of high-energy photons to the photoheating rate is reduced 
due to the weighting by the absorption cross-section am{i^), 
which is a strongly decreasing function of the photon energy. 
Observe that the temperatures (like the neutral fractions) 
obtained in the grey thin simulation asymptote to those ob- 
tained in the mult i- frequency simulation at small distances 
to the ionising source, while the temperatures in the grey 
thick simulation are too high. 

We summarise our discussion of the differences between 
the grey and multi-frequency simulations for the present 
problem by noting that the use of the grey approximation 
leads to neutral fractions and temperatures that generally 
are very different from those obtained in detailed multi- 
frequency simulations. At large optical depths, the neutral 
fractions are systematically too high and the temperatures 
too low due to the lack of spectral hardening. The grey treat- 
ment yields neutral fractions and temperatures that asymp- 
tote to those obtained in the corresponding multi-frequency 
simulation at small distances to the ionising source when 
photoheating rates are computed in the optically thin limit, 
i.e. using Eq. 16. When computing photoheating rates in the 
optically thick limit, i.e. using Eq. 17, the neutral fractions 
and temperatures do not asymptote to the correct values at 
small distances to the ionising source, i.e. the values in the 
multi-frequency simulation. Consequently, when one invokes 
the grey approximation to compute the thermal structure 
of ionised regions, one should compute photoheating rates 
in the optically thin limit. Photoheating rates in the opti- 
cally thick limit should only be employed when considering 
the thermal balance of an ionised region as a whole. Ideally, 
one would perform detailed multi-frequency simulations and 
simply dispense with the grey approximation. 

We now discuss the results of our simulations with 
ttId with respect to those obtained with C^-RAY, CRASH 
and FTTE for the same test problem (Iliev et al. 2006a). 
We note that the simulation with CRASH employed multi- 
ple frequency bins, while the one with ftte was done us- 
ing a single frequency bin and computing photoionisation 
and optically thick photoheating in the grey approximation 
(Alexei Razoumov, private communication). Finally, C^-RAY 
used a hybrid method (Garrelt Mellema, private communi- 
cation): the absorption of ionising radiation was computed 
as a function of frequency, but each photoionisation injected 
the same amount of energy, regardless of the frequency of 
the absorbed photon. This method thus accounts fully for 
the spectral hardening of the radiation but ignores it when 
computing photoheating rates. 

There are noticeable differences in the results obtained 
with these three codes. At large distances from the ionising 
source, i.e. close to and beyond the ionisation front, most 
of these differences may be attributed to differences in the 
multi-frequency implementation, leading to differences in 
the spectral hardening of the emitted blackbody spectrum. 
At these distances, the neutral fractions obtained in our grey 
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simulations agree closely with those obtained with ftte, 
while the neutral fractions obtained in our multi-frequency 
simulations closely agree with (and are, in fact, nearly iden- 
tical to) those obtained with C^-RAY, as expected from our 
discussion above. We note that the fact that the neutral frac- 
tions obtained with CRASH are systematically too large may 
indicate that the radiation field was too poorly sampled (see 
Maselli, Ferrara, & Ciardi 2003, in particular their Fig. 2, 
for a thorough discussion). 

The results, however, show also significant differences in 
the neutral fractions and temperatures close to the ionising 
source, where the gas is close to optically thin and the emit- 
ted blackbody radiation spectrum is not severely deformed 
due to spectral hardening. Some of these differences can be 
attributed to the fact that the different codes employ differ- 
ent expressions for cross-sections, recombination and cooling 
rates. As demonstrated in Iliev et al. (2006a) (their Fig. 4), 
different recombination and cooling rates may only account 
for differences in the neutral fraction and temperature of at 
most < 10%. We have verified this by employing the rates 
used with the different codes (Table 2 in Iliev et al. 2006a) 
in simulations with ttId. 

Most of the differences close to the ionising source may 
instead be traced back to the use of different assumptions 
underlying the computation of the photoheating rates. In 
fact, the temperatures obtained with CRASH are in very good 
agreement with the temperatures obtained in our multi- 
frequency and grey thin simulations, while the temperatures 
obtained by ftte and C^-RAY are in excellent agreement 
with the temperatures in our grey thick simulation. 



5.2.3 HII region in gas containing hydrogen and helium 
with an evolving temperature 

Finally, we test the ability of ttId to accurately compute 
the ionisation and temperature structure in gas containing 
both hydrogen and helium by comparing results obtained 
with ttId with results obtained with the photoionisation 
code CLOUDY (version 08.00; last described by Ferland et al. 
1998) for the same test problem. Note that cloudy assumes 
ionisation equilibrium. 

As before we consider an ionising source with black- 
body spectrum of temperature T = 10^ K with an ionising 
luminosity of 5 x 10^^ s~^ in gas of density nu = 10~^, 
but we now assume X — 0.75 and Y = 1 — X . cloudy in- 
cludes considerably more physics than ttId. To facilitate a 
direct comparison, we therefore keep the setup of the sim- 
ulations as simple as possible: we assume that there is no 
radiative coupling between hydrogen and helium, i.e., pho- 
tons emitted due to recombination of helium do not lead to 
ionisations of hydrogen, and compute recombinations in the 
Case A (Sec. 3.1.2) limit. 

Fig. 5 shows the ionised fractions and temperatures 
computed with ttId at time t — 2000 Myr. It also shows the 
results obtained with cloudy, which correspond to a time 



^'^ Maselli, Ciardi, & Kanekar (2009) have repeated this test with 
a more recent version of CRASH with improved sampling of the 
Monte Carlo photon field. They find slightly larger temperatures 
(their Fig. 3), which further improves the agreement with the 
temperatures found with ttId. 
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Figure 5. Comparison of the results from a simulation with ttId 
(solid curves) with those from a simulation with CLOUDY (dashed 
curves) for the HII region test problem described in Sec. 5.3 (as- 
suming X = 0.75, y = 1 — X). Note that the results computed 
with CLOUDY correspond to equilibrium {t — oo), while the re- 
sults computed with ttId correspond to a time t = 2000 Myr 
(which is much larger than the recombination time on which the 
gas reaches thermal and ionisation equilibrium). The results ob- 
tained with ttId are nearly indistinguishable from those obtained 
with cloudy for all radii at which equilibrium has been reached. 



t ^ oo. The agreement between ttId and cloudy is excel- 
lent. The results only differ at radii where (in the simula- 
tion with ttId) equilibrium has not yet been reached. Thus, 
ttId yields equilibrium ionisation and temperature profiles 
that are almost indistinguishable from those obtained with 
a well-tested and widely employed state-of-the-art photoion- 
isation code. 



5.3 Test 2: HII region expansion. TRAPHIC 

In this section we apply our new implementation of traphic 
to compute the evolution of the ionisation state and tem- 
perature around an ionising source surrounded by gas of 
constant density. This idealised test problem captures the 
main characteristics of a thermally coupled RT simulation 
that we wish to verify: conservation of the number of ion- 
ising photons, which ensures that the final ionised region 
attains the correct size, and conservation of the associated 
energy, which, together with an accurate implementation of 
the relevant cooling processes, ensures that the ionised re- 
gion settles into the correct thermal structure. The physical 
parameters for the test are identical to those employed to 
obtain the reference solutions presented in Sec. 5.2.2 but, 
for definiteness, we repeat the problem description here. 

We consider an ionising source embedded in a homo- 
geneous hydrogen-only density field with number density 
riH — 10~^ cm~^. The source emits = f du JV^^v) — 
5 X 10^^ ionising photons s~^ with a blackbody spectral 
shape J\fj{iy) corresponding to a blackbody temperature 
Tbb = 10^ K. The test described here is identical to Test 1 
in Paper I, except that now the gas temperature is allowed 
to vary due to heating and cooling processes as described 
in Sec. 3.2 (with Compton cooling off the redshift z = 
cosmic microwave background included) and that collisional 
ionisation is included. The gas is assumed to have an initial 
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hydrogen ionised fraction t/hii = 1.2 x 10~^ (approximately 
corresponding to the ionised fraction imphed by cohisional 
ionisation equihbrium at the temperature T = 10^ K). Its 
initial temperature is set to 10^ K. 

First, in Sec. 5.3.1, we consider the case in which radi- 
ation is transported using a single frequency bin in the grey 
optically thin approximation in pure hydrogen gas. We em- 
ploy the grey approximation to allow for a more direct com- 
parison with the results presented in Paper L In Sec. 5.3.2 
we will also briefly discuss the performance of TRAPHIC in 
multi-frequency simulations in gas containing both hydrogen 
and helium. 



5.3.1 HII region expansion: grey thin, hydrogen-only 

The numerical realisation of the initial conditions is simi- 
lar to that used for Test 1 in Paper I. The ionising source 
is located at the centre of a simulation box with side 
length Lbox = 13.2 kpc. The box boundaries are photon- 
transmissive, i.e., photons leaving the box are lost from the 
computational domain. We assign each SPH particle a mass 
m — nnmnLi^o^/^spUj where A^sph is the total number 
of SPH particles. The positions of the SPH particles are 
chosen to be glass- like (e.g.. White 1996). Glass- like initial 
conditions imply a more regular distribution of particles in 
space when compared to that obtained from a Monte Carlo 
sampling of the density field. The SPH smoothing kernel is 
computed and the SPH densities are found using the SPH 
formalism implemented in GAD get- 2, with A^ngb = 48. 

Photons are transported using a single frequency bin 
assuming the grey approximation in the optically thin 
limit. We therefore employ a photoionisation cross-section 
(cr7Hi) = 1.63 X 10~^^ cm^ (Sec. 3.1) and assume that 
each photoionisation adds (chi) = 6.32 eV to the thermal 
energy of the gas (Sec. 3.2). The RT time step is set to 
Atr = 10~^ Myr to facilitate a comparison to Test 1 in Pa- 
per I. For the same reason, we limit ourselves to solving the 
time-independent RT equation and propagate photons dur- 
ing each time step only from a given particle to its direct 
neighbours. All simulations presented in this section employ 
NspK — 64^ SPH particles, which are evolved for a total of 
500 Myr. Some of our simulations employ the resampling 
technique introduced in Paper I to reduce artifacts due to 
the particular setup of the initial conditions. Briefly, each 
SPH particle is, within its spatial resolution element whose 
size is determined by the diameter of the SPH kernel, 2/i, 
from time to time^^ offset randomly from its initial position. 
For comparison, we repeat all simulations without employ- 



The particle distribution is resampled every 10th RT time step. 
Our results are insensitive to the precise frequency with which 
the resampling is applied. We note that the choice for the resam- 
pling frequency is problem-dependent and hence the resampling 
frequency must usually be determined experimentally using con- 
vergence tests. In simulations with many sources, in which SPH 
particles receive photons from many different directions, artefacts 
due to the particular arrangement of SPH particles are typically 
much less prominent, as discussed in Paper I (Sec. 5.3.3; see also, 
e.g., the related discussion on cell randomization in Trac & Cen 
2007). Hence, realistic simulations will typically not require re- 
sampling. 



ing this technique. We perform simulations with different 
angular resolutions. Figs. 6 and 7 show our results. 

In Fig. 6 we present slices through the centre of the 
simulation box showing the neutral fraction (top row) and 
temperature (bottom row) at time^^ t = 100 Myr. In each 
row, the three left-most panels show results from simulations 
with angular resolution A^c = 8, 32 and 128 and no resam- 
pling of the particle positions applied and the right-most 
panel shows results from a simulation with angular resolu- 
tion Ac = 32 and resampling of the particle positions every 
10th RT time step. In each panel we indicate, as a point of 
reference, the analytical approximation for the position of 
the ionisation front (Eq. 31) by a dash-dotted circle. 

Interior to the ionisation front the gas is highly ionised 
and photo-heated to typical temperatures T 1.5 x 10^ K 
(with maximum temperatures T ^ 2 x 10^ K). The runs 
that did not employ the resampling show slight deviations 
from the expected spherical shape which depend on the an- 
gular resolution. As discussed in Paper I, the deviations 
are caused by the particular arrangement of the SPH par- 
ticles. Reducing this particle noise, which is strongest when 
A^c ~ A^ngb, was the motivation for introducing the resam- 
pling technique. Indeed, the distribution of neutral fractions 
and temperatures from the simulation that employed the re- 
sampling of the density field is spherically symmetric to a 
high degree. 

In Fig. 7 we compare the median profiles of the neu- 
tral fraction (left-hand panel) and the temperature (right- 
hand panel) at time t = 100 Myr obtained from the three- 
dimensional simulations with traphic (solid curves with er- 
ror bars) to the reference simulation obtained with our 1-d 
RT code ttId (dashed curves). The results of all simula- 
tions are in excellent agreement with the reference result. 
The small deviations that are present very close to the ion- 
ising source and in regions where the profile gradients are 
steep are due to the finite spatial resolution (indicated with 
horizontal error bars in the top left corners of the panels). 
The effect of resampling in reducing noise can most clearly 
be seen when comparing the simulations with angular resolu- 
tion A^c = A^ngb — 32 with each other (middle panels). Note 
that for the simulation with the highest angular resolution 
that we have considered here (Ac = 128), the resampling 
slightly reduces the agreement with the reference simula- 
tion because it introduces additional scatter. This scatter is 
consistent with the spatial resolution employed. 



5. 3. 2 HII region expansion: multi- frequency, hydrogen and 
helium 

Next we demonstrate the ability of TRAPHIC to accurately 
solve the present multi-frequency problem in gas of primor- 
dial composition (i.e., in the presence of helium) and using 
multiple frequency bins. For brevity, we discuss only a single 
simulation with particle number A" = 64^ and angular reso- 



The reason why we do not show the slices at the end of the 
simulations, i.e. at time t = 500 Myr, as we did in the correspond- 
ing Test 1 in Paper I, is that the simulation box is slightly too 
small to contain the whole ionised sphere at this time (because of 
the smaller photoionisation cross-section that is employed here). 
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Figure 6. Test 2 (using traphic). Neutral fraction (top row) and temperature (bottom row) at time t = 100 Myr in a slice 
through the centre of the simulation box. From left to right: angular resolution A^c = 8, 32, 128 (all without resampling) and 32 
(with resampling of the particle positions after every 10th RT time step, as indicated by the letter 'R' in the panel titles). The 
dot-dashed circles indicate the position of the ionisation front, calculated using the analytical approximation (Eq. 31). Contours 
show neutral fractions of ?7hi = 0.9, 0.5, log^^Q ?7hi = —1, —1.5, —2, —2.5, —3, —3.5 and temperatures log]^Q(T/K) = (3, 4, 4.2) (from 
the outside in). The crosses in the top row indicate the spatial resolution (2h). 
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Figure 7. Test 2 (using traphic). Scatter plots and profiles of the neutral (ionised) fraction and temperature at time 
t = 100 Myr for simulations with angular resolution Nc = 8 (left), 32 (middle) and 128 (right). Top row: No resampling. 
Bottom row: Resampling of the particle positions after every 10th RT time step (indicated by the letter 'R' in the legends). 
Each dot represents the neutral fraction (ionised fraction, temperature) of a single SPH particle (only a randomly chosen 
subset of 10% of all particles is shown). Solid curves show the median neutral fraction (ionised fraction, temperature) in 
spherical bins around the ionising source. The vertical error bars enclose 68.3% of the particles in each bin. Dashed curves 
indicate the reference solution obtained with our 1-d RT code ttId. The horizontal error bars in the upper left corners 
indicate the spatial resolution. The results of all simulations are in excellent agreement with the reference solution. Without 
the resampling, the results are noisier if A^c ~ ^ngb (top middle panel). 



© RAS, MNRAS 000, 1-23 



TRAP HI C - multi-frequency and thermal coupling 17 
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Figure 8. Test 2 (using traphic). Multi-frequency transport and 
effects of the inclusion of helium (Sec. 5.3.2). Neutral hydrogen 
fraction ryni {top panels) and temperature (bottom panels) at time 
t = 100 Myr in a slice through the centre of the simulation box. 
The left and right panels show, respectively, results from simu- 
lations of an ionising source in gas of pure hydrogen (X = 1) 
and gas with primordial abundances (X = 0.75, Y = 1 — X). 
Except for this difference in the helium abundance, both simu- 
lations are identical; in particular, they both make use of five 
frequency bins to transport the emitted 10^ K blackbody pho- 
tons. The panels can be compared with the left-most panels of 
Fig. 6, which show results from a simulation that is identical ex- 
cept that it used a single frequency bin and the grey optically 
thin approximation. As in Fig. 6, the dot-dashed circle indicates 
the position of the ionisation front, calculated using the analyt- 
ical approximation (Eq. 31). Contours show neutral fractions of 
^Hl = O-Q, 0.5 , logj^Q ?7hi = —1, —1.5, —2, —2.5, —3, —3.5 and tem- 
peratures logio(r/K) = (3,4,4.2,4.4) (from the outside in). The 
colour coding is the same as in Fig. 6. The crosses indicate the 
spatial resolution {2h). 



lution Nc = 8. We have verified that simulations with other 
choices for these parameters show the expected behaviour. 

We perform two simulations. The first simulation as- 
sumes a hydrogen mass fraction X = 1. The second simu- 
lation assumes a hydrogen mass fraction X = 0.75 and a 
helium mass fraction Y = 1 — X . We set the initial ionised 
helium fractions to zero, ?7Heii = rjueiii — 0. All other physi- 
cal parameters are as in the previous section. For both sim- 
ulations we use the same number of frequency bins, Ni, — 5 
(starting at 13.6 eV, 24.6 eV, 35.5 eV, 54.4 eV and 75.0 eV, 
with the last bin extending to infinity) . The photoionisation 
cross-section and excess energy associated with each bin are 
obtained from averaging over a blackbody spectrum of tem- 
perature 10^ K, assuming the optically thin limit (Eqs. 9 
and 16). 

The motivation behind our choice to use a small num- 
ber of frequency bins is that in realistic simulations that 



will be computationally more expensive, limited resources 
will require the usage of as few frequency bins as possible. 
Our results below show that a number as low as five (and 
perhaps even as low as three, see Sec. 5.4) may be sufficient 
to capture the main effects associated with multi-frequency 
radiation transport. 

Fig. 8 shows the neutral hydrogen fractions (top) and 
temperatures (bottom) in a slice through the centre of the 
simulation box for the simulation without (left) and with 
(right) helium. The panels can be compared with the left- 
most panels of Fig. 6, which show results from an identical 
simulation except that it used a single frequency bin and 
the grey, optically thin approximation. The effect of spectral 
hardening is most visible in the panels showing the temper- 
ature, with the multi-frequency simulations showing a sub- 
stantial pre- heating ahead of the ionisation front. Interest- 
ingly, the simulation that includes helium shows slightly less 
pre-heating (and pre-ionisation) than the one that assumes 
pure hydrogen. 

The solid curves in Fig. 9 show median profiles of the 
species fractions and temperatures for both simulations. For 
comparison, the (converged) reference results obtained with 
ttId are shown by dashed curves. The results of the sim- 
ulations with TRAPHIC are in excellent agreement with the 
reference result, both with and without the inclusion of he- 
lium. The small deviations that are present very close to the 
ionising source and in regions where the profile gradients are 
steep are due to the finite spatial resolution (indicated with 
horizontal error bars). The fact that ?7Heii shows reduced 
scatter is probably because its value is not free but depends 
on ?7Hei and ?7Heiii according to Eq. 5. 

5.4 Test 3: Cosmological reionisation 

In this section we use our thermally coupled implemen- 
tation of TRAPHIC to repeat Test 4 of the cosmological 
RT code comparison project (Iliev et al. 2006a) that we 
have discussed in Paper I for the case of fixed temperature 
(T = 10^ K). This test involves the simulation of the evo- 
lution of ionised regions around multiple sources in a static 
cosmological density field at redshift z ~ 8.85 and was de- 
signed to resemble important aspects of state-of-the-art sim- 
ulations of the epoch of reionisation. In contrast to our Test 
4 simulations in Paper I, we will here compute the evolution 
of the temperature along with that of the ionisation state 
of the gas. We will perform both simulations assuming the 
grey approximations (optically thin and optically thick) and 
using multiple frequency bins and we will discuss the origin 
of the differences in the results that these simulations yield. 

The setup of this test is identical to that of Test 4 in 
Paper I, to which we refer the reader for a detailed de- 
scription. Briefly, the initial conditions are provided by a 
snapshot (at redshift z ^ 8.85) from a cosmological N-body 



15 While this statement is certainly true for the present test prob- 
lem, we caution that the answer to the question of how many fre- 
quency bins are sufficient will be problem-dependent. Hence, the 
minimum number of frequency bins that can be employed while 
still capturing the main physical effects must be determined by 
performing explicit convergence studies for the particular problem 
at hand (see also, e.g., McQuinn et al. 2009). 
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Figure 9. Test 2 (using traphic). Multi-frequency transport and effects of the inclusion of helium (Sec 5.3.2). Scatter plots and profiles 
of species fractions and temperature at time t = 100 Myr. The left and right panels show, respectively, results from simulations of an 
ionising source in gas of pure hydrogen (X = 1) and gas with primordial abundances (X = 0.75, Y = 1 — X). Except for this difference in 
the helium abundance, both simulations are identical; in particular, they both make use of five frequency bins to transport the emitted 
10^ K blackbody photons. The panels could be compared with the left-most panels of Fig. 7, which show results from a simulation that 
is identical except that it used a single frequency bin and the grey, optically thin approximation. As in Fig. 7, each dot represents the 
species fraction (grey: ?7hi, red: ?7hii, blue: ?7Hel5 green: ?7Hell5 yellow: ?7Helll) or temperature of a single SPH particle (only a randomly 
chosen subset of 10% of all particles is shown). Solid curves show the median neutral fraction (ionised fraction, temperature) in spherical 
bins around the ionising source. The vertical error bars enclose 68.3% of the particles in each bin. Dashed curves indicate the reference 
solution obtained with our 1-d RT code ttId. The horizontal error bars in the upper left corners indicate the spatial resolution {2h). 



and gas- dynamical uniform- mesh simulation. The simula- 
tion box is Lbox = 0.5 comoving Mpc on a side, uni- 
formly divided into A^ceii = 128^ cells. We Monte Carlo 
sample this input density field to replace the mesh cells with 
Nspu = A^ceii = 128^ SPH particles. The gas is assumed to 
be initially fully neutral and to have an initial temperature 
of T = 10^ K. The ionising sources are chosen to corre- 
spond to the 16 most massive halos in the box. They are 
assumed to have blackbody spectra Bu{i^,Thh) with tem- 
perature Tbb = 10^ K. The ionising photon production rate 
is taken to be constant and all sources are switched on at the 
same time. The box boundaries are photon-transmissive. 

We perform three RT simulations to solve the time- 
independent RT equation, all with an angular resolution of 
A^c = 32 (and setting A^ngb = 32). We have demonstrated 
in Paper I (Test 4) that for the current problem this an- 
gular resolution is sufficiently high to obtain converged re- 
sults. To facilitate a direct comparison with the correspond- 
ing simulation in Paper I, we employ the same time step 
Atr = 10~^ Myr and transport photons only over a single 
inter-particle distance per time step. We note that the cur- 
rent simulations do not employ the resampling technique to 
suppress noise in the neutral fraction caused by the par- 
ticular realisation of the SPH density field. As discussed in 
Paper I, for the present test this noise is small. For definite- 
ness we mention that all simulations include collisional ioni- 
sation and all relevant cooling processes (including Compton 



cooling off the z = 8.85 cosmic microwave background), em- 
ploying the rates listed in Table 1. 

In the first two of the three simulations we transport 
radiation using a single frequency bin, employing the grey 
photoionisation cross-section {a^m) = 1.63 x 10~^^ cm^ 
(Sec. 3.1). The difference between these two simulations is 
in the computation of the photoheating rates used to evolve 
the gas temperatures. For one simulation (traphic thin) we 
compute photoheating in the optically thin limit, assuming 
that each photoionisation adds (chi) = 6.32 eV to the ther- 
mal energy of the gas (Sec. 3.2). In the other simulation 
(traphic thick) we compute photoheating in the optically 
thick limit, assuming that each photoionisation on average 
adds (em''^) = 16.01 eV to the thermal energy of the gas 
(Sec. 3.2). 

The third simulation (traphic) differs from the first 
two in that we transport photons using A^^^ = 3 frequency 
bins (starting at 13.6 eV, 35 eV and 50 eV, with the last 
bin extending to infinity). The photoionisation cross-section 
and the excess energy associated with each bin are obtained 
from averaging over a blackbody spectrum of temperature 
10^ K, assuming the optically thin limit (Eqs. 9 and 16). As 
in the previous section, our choice in favour of a very small 
number of frequency bins has purely practical reasons: com- 
putational efficiency. While we could have performed this 
relatively small test simulation at higher spectral resolution, 
we anticipate that applications of traphic to large simula- 
tions of reionisation will generally require us to choose a 
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Figure 10. Test 3. Neutral fraction in a slice through z = Lbox/2 at time t = 0.2 Myr. From left to right: TRAPHIC thin (assuming grey 
optically thin photoheating rates), TRAPHIC thick (assuming grey optically thick photoheating rates), TRAPHIC (using three frequency 
bins), C^-RAY, CRASH and FTTE. Contours show neutral hydrogen fractions ?7hi = 0.9,0.5, log?7Hi = —1, —3 and —5, from the outside in. 
The results obtained with TRAPHIC thick are in excellent agreement with those obtained with FTTE. They are also in excellent agreement 
with the results obtained with C-^-RAY in highly ionised regions, where the neutral fraction is unaffected by spectral hardening. The 
small differences in the neutral fractions obtained with TRAPHIC thick, TRAPHIC thin and TRAPHIC are mostly due to differences in the 
recombination rate, caused by differences in the gas temperatures (see Fig. 11). 
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Figure 11. Test 3. Temperature in a slice through z = Lbox/2 at time t = 0.2 Myr. From left to right: TRAPHIC thin (assuming 
optically thin photoheating rates), TRAPHIC thick (assuming optically thick photoheating rates), TRAPHIC (using three frequency bins), 
C^-RAY, CRASH and FTTE. Contours show temperatures \ogiQ(T /K) = 3, 4, 4.2, 4.4 and 4.6, from the outside in. Most of the morphological 
differences may be attributed to differences in the spectral hardening of the ionising radiation, with the multi-frequency codes TRAPHIC, 
C^-RAY and CRASH yielding a substantial amount of pre-heating and the monochromatic (grey) codes TRAPHIC thin, TRAPHIC thick and 
FTTE yielding sharp transitions between the hot ionised and the cold neutral phases. The differences in the maximum gas temperatures 
are mainly due to photoheating being computed in the optically thick limit (traphic thick, C^-RAY, ftte), the optically thin limit 
(traphic thin) or using multiple frequency bins (traphic, crash). 



number of frequency bins as small as possible. Using a small 
number of frequency bins in the present test should thus 
give results that more closely resemble the results of future, 
larger simulations. 

Figs. 10-12 show our results. Fig. 10 shows images of the 
neutral fraction in slices through the centre of the simula- 
tion box at time t = 0.2 Myr (our conclusions also hold for 
other times). The individual panels show results obtained 
with TRAPHIC thin, traphic thick and traphic. For ref- 
erence, we also show the results obtained with other RT 
codes for the same test problem as published in Iliev et al. 
(2006a). Neutral fraction contours are shown to facilitate the 
comparison. While the simulation with CRASH treated the 
present problem by performing a multi-frequency computa- 
tion, the simulation with ftte, as our simulation traphic 
thick, solved it in the grey approximation using optically 
thick photoheating rates. Finally, C^-RAY employed a hybrid 
method that treats the transport of radiation with multiple 
frequency bins but computes photoheating rates in the grey 
(optically thick) approximation. 



The differences in the neutral fractions are generally 
small. The simulations that employ photoheating rates in 
the optically thick limit (ftte, c^-ray, traphic thick) 
yield smaller minimum neutral fractions than the simu- 
lations that compute photoheating rates in the optically 
thin limit (traphic thin) or using multiple frequency bins 
(crash, traphic). This is the result of lower recombina- 
tion rates caused by the higher temperatures these simula- 
tions yield^^ (see Fig. 11). The regions with low ionisation 
(^Hi > 0.5) found with traphic are slightly smaller than 
those found with C^-RAY and CRASH, which indicates that 
three frequency bins are not sufficient for obtaining highly 
accurate mult i- frequency solution. Still, the simulations with 
TRAPHIC seem to capture the main effects (see the discussion 
on pre-heating below). 

Fig. 11 shows images of the gas temperature in slices 

As noted earlier, the main reason why CRASH finds significantly 
larger neutral fractions may be an insufficient sampling of the 
photon field, see, e.g.. Fig. 2 in Maselli, Ferrara, Sz Ciardi 2003. 
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Figure 12. Test 3. Histograms of the temperature at time t = 
0.2 Myr from simulations with different RT codes. At low temper- 
atures the differences in the shapes of the histograms are mainly 
due to differences in the code-specific treatment of spectral hard- 
ening. The differences exhibited at high temperatures are mainly 
due to photoheating being computed in the optically thick limit 
(traphic thick, C^-ray, ftte), the optically thin limit (traphic 
thin) or using multiple frequency bins (CRASH, traphic). 

through the centre of the simulation box that correspond to 
the images of the neutral fraction shown in Fig. 10. There 
are significant differences in both the morphologies of the 
photo-heated regions and the typical temperatures attained 
by the photoionised gas between the different simulations. 
Outside the ionisation fronts, these differences can mostly 
be attributed to differences in the spectral hardening of the 
ionising radiation. C^-RAY, CRASH and traphic all yield a 
substantial pre-heating of the gas ahead of the ionisation 
fronts. This pre-heating is not seen in the simulations with 
TRAPHIC thin, TRAPHIC thick and ftte since they all as- 
sume the grey approximation. In Sec. 5.2.2 we have already 
discussed, for the same set of codes, the differences between 
a multi-frequency treatment and its grey approximations in 
idealised simulations of the evolution of a single, spherically 
symmetric, ionised region. The results here are in close qual- 
itative agreement with that discussion. 

The results obtained with the different codes also ex- 
hibit significant variations in the gas temperature in regions 
well inside the ionisation fronts. While CRASH, traphic and 
TRAPHIC thin yield typical temperatures^^ of T 2 — 3 x 
10^ K, the typical temperatures obtained with C^-RAY, ftte 
and TRAPHIC thick are T ^ 6 x 10^ K, i.e., substantially 
higher. Note that there is also disagreement between codes 
which incorporate the detailed treatment of multi-frequency 



Maselli, Ciardi, &; Kanekar 2009 have repeated this test with 
a more recent version of CRASH with improved sampling of the 
Monte Carlo photon field. They find slightly larger (20 - 30% at 
t = 0.05 Myr) peak temperatures (their Fig. 5), which improves 
the agreement with the peak temperatures found with traphic. 



radiation (c^-RAY, CRASH, traphic) and between codes in 
which the radiation is treated in the grey approximation 
(ftte, traphic thin, traphic thick). Spectral hardening 
may thus only provide an explanation for part of the differ- 
ences between the simulations. 

We recall that in Sec. 5.2.2, where we simulated the 
evolution of a single, spherically symmetric, photoionised 
region, we found qualitatively similar differences between 
the results obtained with C^-RAY, CRASH and ftte. In the 
(nearly optically thin) region close to the ionising source, 
the simulations that employed C^-RAY and ftte yielded gas 
temperatures that were substantially larger than those in 
the simulation that employed CRASH. By comparing with 
results obtained with our 1-d RT code ttId, we were able 
to explain most of these temperature differences in terms of 
differences in the assumptions underlying the computation 
of photoheating rates. The results presented in Fig. 11 are 
another manifestation of this explanation. 

In Fig. 12 we compare histograms of the temperature at 
time t = 0.2 Myr. For the simulations with C^-RAY, CRASH 
and FTTE we have computed these histograms directly from 
the A^ceii = 128^ values of temperatures published in Iliev 
et al. (2006a). For the simulations with traphic, traphic 
thin and TRAPHIC thick we assigned the temperatures to a 
corresponding uniform mesh with A^ceii = 128^ cells using 
SPH interpolation before we computed the histograms. 

The histograms provide a quantitative confirmation 
of our qualitative discussion above. The simulations with 
TRAPHIC and TRAPHIC thin yield, in close agreement with 
the simulations performed with CRASH, typical temperatures 
ofT?^2 — 3xl0^K. On the other hand, the simulations per- 
formed with C^-RAY, FTTE and traphic thick closely agree 
on typical temperatures of T 6 x 10^ K. The differences 
between the histograms at low values of the temperature 
are mostly caused by the differences in spectral hardening. 
Due to the pre-heating of gas ahead of the ionisation fronts 
obtained with the multi- frequency codes TRAPHIC, C^-RAY 
and CRASH, the number of cells that are still at their ini- 
tial temperature T = 100 K is much smaller than found by 
FTTE, TRAPHIC thin and traphic thick, which only employ 
a single frequency bin. 

In summary, in this section we have repeated the simu- 
lation of the expansion of multiple ionised regions in a cos- 
mological density field that we discussed in Test 4 in Paper I, 
but this time we explicitly computed the evolution of the 
temperature of the photoionised gas. We performed three 
simulations: one simulation computed photoheating in the 
grey, optically thin limit (traphic thin) , one computed pho- 
toheating in the grey, optically thick limit (traphic thick) 
and one employed multiple (i.e., three) frequency bins. All 
three simulations showed only small differences in the neu- 
tral fractions when compared with each other, which could 
be plausibly explained with the differences in the temper- 
atures. The temperatures differed due to the computation 
of photoheating rates in different limits (grey optically thin, 
grey optically thick and multi- frequency). In particular, em- 
ploying the grey approximation in the optically thick limit 
yields temperatures that are typically too large by factors 
2 — 3. We also compared the results of our thermally coupled 



See footnote 17 
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simulations with results obtained with other RT codes for 
the same test problem (Iliev et al. 2006a). We found very 
good agreement between these and our results when com- 
paring simulations that employed similar assumptions for 
computing photoionisation and photoheating rates. 

The results of the multi-frequency simulation (neutral 
fractions and temperatures interpolated to a 128^ uniform 
mesh using SPH interpolation) are available for download 
at the website of the Cosmological RT Code Comparison 
Project (Iliev et al. 2006a; Iliev et al. 2009). 



6 SUMMARY 

In this work we described an extension of the implemen- 
tation of TRAPHIC, the radiative transfer (RT) method for 
use with Smoothed Particle Hydrodynamics (SPH) simula- 
tions that we have introduced in Pawlik & Schaye (2008, 
Paper I), in a modified version of the SPH code GADGET 
(Springel 2005). The new implementation of traphic can 
be used to solve multi-frequency RT problems in primor- 
dial gas consisting of hydrogen and helium. It also allows 
for the computation of the non-equilibrium evolution of the 
gas temperature due to photoheating and radiative cooling. 

As part of the new implementation we introduced a 
numerical method that allows us to accurately compute the 
coupled evolution of the ionisation balance and temperature 
of gas parcels exposed to ionising radiation and that works 
independently of the size of the chosen RT time step. This 
decoupling of the RT time step from the time scales that gov- 
ern the evolution of the species fractions and temperatures, 
i.e. the decoupling from the ionisation, recombination and 
radiative cooling time scales, is an important pre-requisite 
for performing efficient RT simulations. The alternative, a 
RT time step limited by the values for the ionisation, re- 
combination or cooling time scales, could quickly become 
computationally infeasible since these time scales may be- 
come very small. 

We discussed the performance of the new, thermally 
coupled implementation of traphic in three-dimensional 
multi-frequency RT test simulations of spherically symmet- 
ric expanding HII regions and non-spherical expanding HII 
regions in a cosmological reionisation setting. We treated 
the multi-frequency radiation both using a single frequency 
bin in the grey approximation and using multiple frequency 
bins. We distinguished two types of grey approximations by 
computing photoheating both in the optically thin and op- 
tically thick limits. We compared the results of our test sim- 
ulations to results obtained with other RT codes for iden- 
tical test problems. We found excellent agreement in the 
morphologies and gas temperatures of the photoionised and 
photoheated regions when comparing simulations that em- 
ployed similar assumptions for computing photoionisation 
and photoheating rates. 

We used the new implementation to demonstrate and 
pinpoint the differences in the results obtained from grey 
simulations and simulations that use multiple frequency 
bins. Close to and ahead of ionisation fronts these differ- 
ences are mostly due to the spectral hardening of the ra- 
diation field caused by the dependence of the absorption 
cross-section on photon energy. Spectral hardening signifi- 
cantly increases the widths of ionisation fronts and implies a 



substantial preheating of the gas ahead of them. Additional 
significant differences between the simulations are caused 
by the choice of the limit in which grey photoheating rates 
are computed. Simulations that use grey photoheating rates 
computed in the optically thick limit yield typical gas tem- 
peratures that are too large by factors 2 — 3 when compared 
with the exact multi-frequency solution. Simulations that 
use grey photoheating rates computed in the optically thin 
limit yield typical gas temperatures that asymptote to the 
multi-frequency result with decreasing distances from the 
ionising sources. 

The additions presented in this work are crucial for ap- 
plications of traphic to simulations of the (re-) ionisation 
of both hydrogen and helium that also wish to account for 
the preheating of gas ahead of ionisation fronts due to spec- 
tral hardening. We plan to perform such simulations in the 
future. Note that we have limited our considerations to RT 
simulations on pre-computed static density fields. But pho- 
toheating increases also the gas pressure and hence affects 
the hydrodynamical evolution of the gas. An important goal 
for the future will therefore be to present an implemen- 
tation of TRAPHIC that allows one to perform radiation- 
hydrodynamical simulations. 
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APPENDIX A: A NEW TREATMENT OF 
ABSORPTIONS BY VIRTUAL PARTICLES 

In this appendix we show that the treatment of virtual par- 
ticles (ViPs) in the implementation of traphic that we have 
used to perform the (hydrogen-only) simulations published 
in Pawlik & Schaye (2008, hereafter Paper I) and that we 
will refer to as the old implementation, results in a tempo- 
rary underestimate of the neutral hydrogen fraction just be- 
hind evolving ionisation fronts in simulations that use a high 
angular resolution. We will show that this underestimate is 
absent in simulations that employ our new implementation 
(presented in the current work). Moreover, in simulations 
that employ this new implementation, the numerical scat- 
ter in the neutral hydrogen fraction is significantly reduced. 
For clarity of the presentation and because we will compare 
results obtained with the new implementation with results 
obtained with the old implementation presented in Paper I 
that lacked the treatment of helium, we will assume that 
photons are transported in gas that consists only of hydro- 
gen (i.e., X = 1). Our discussion generalizes straightfor- 
wardly to the transport of photons in gas consisting of both 
hydrogen and helium. 

The number of hydrogen-ionising photons a ViP ab- 
sorbs depends on its neutral hydrogen density. As explained 
in Paper I, the computation of this number is performed in 
exactly the same manner as for SPH particles. The only dif- 
ference between the treatment of photons absorbed by SPH 
particles and ViPs is that the latter distribute the absorbed 
photons amongst their SPH neighbours. For this distribu- 
tion of absorbed photons one must specify the fraction of 
the total that is given to each of the SPH neighbours. In 
the old implementation of traphic this fraction was taken 
to be proportional to the value of the SPH kernel W of the 
distributing ViP at the position of the SPH neighbour. In 
the new version this fraction is taken to be proportional to 
the contribution of the SPH neighbour to the SPH estimate 
of the ViP's neutral hydrogen density. 

The old treatment of ViPs results in an underestimate 
of the simulated non-equilibrium neutral hydrogen fractions. 
Fig. Al serves to demonstrate this. Its panels show the neu- 
tral and ionised hydrogen fractions around a single ionising 
source in a homogeneous hydrogen-only medium at times 
t = 30, 100 and 500 Myr (from left to right) obtained with 
the old (first and third rows) and new (second and fourth 
rows) implementation. The setup and parameters for the 
simulations presented here are identical to the setup and 



© RAS, MNRAS 000, 1-23 



TRAPHIC - multi-frequency and thermal coupling 23 



parameters used for the A^sph = 64^, A^c = 128 simulation 
presented in Test 1 in Sec. 5.3.1 of Paper I. In addition to the 
neutral (grey dots) and ionised (light red dots) fractions of 
each particle, Fig. Al shows the median neutral (black solid 
curves) and ionised (red solid curves) fractions in spheri- 
cal bins, which are compared to the exact solution obtained 
with our 1-d RT code ttId (Sec. 5.3; dashed curves of the 
corresponding colour). The error bars indicate the 68.3% 
confidence intervals in each bin. For each implementation 
we have performed simulations both with and without re- 
sampling the density field, as indicated by the presence or 
absence of the letter 'R' in the panel titles. 

In the simulations employing the old implementation of 
TRAPHIC the neutral hydrogen fractions at times t = 30 and 
100 Myr are underestimated at radii slightly smaller than 
the radius of the ionisation front. In the simulations that 
employ the new implementation this underestimate is no 
longer present, thanks to the new manner in which the pho- 
tons absorbed by ViPs are distributed. At ^ = 500 Myr, i.e. 
when the ionised region has (nearly) reached its equilibrium 
size, the underestimate is also absent in the simulations that 
employ the old implementation. However, at this time these 
old simulations still exhibit an increased scatter around the 
median when compared to the corresponding snapshots from 
the simulations that employ the new implementation. 

The underestimate of the neutral hydrogen fraction just 
behind evolving ionisation fronts in simulations employing 
the old implementation is caused by the fact that in this 
implementation the distribution of the photons absorbed by 
ViPs does not respect the spatial distribution of the neu- 
tral gas in their surroundings. It mainly affects the neu- 
tral fraction of particles close to evolving ionisation fronts, 
because the number of photons absorbed and subsequently 
distributed by ViPs near the ionisation front is significantly 
larger than the number of photons that are absorbed by 
the SPH particles behind the ionisation front and because 
the ViPs distribute the absorbed photons irrespective of the 
neutral hydrogen mass with which the corresponding SPH 
particles contributed to the computation of its neutral hy- 
drogen density. 

We did not notice the described temporary underes- 
timate of the neutral fraction just behind non-equilibrium 
ionisation fronts in the simulations that we have presented 
in our original publication (Paper I), since there we only dis- 
cussed profiles of the neutral fraction at t = 500 Myr. The 
reason why we limited ourselves to discussion of equilibrium 
results in that publication, was that we were still lacking ac- 
curate non-equilibrium reference solutions at that time (our 
1-d reference RT code ttId was still under development). 
The discovery of the underestimate of the neutral fraction 
was triggered by scatter plots of the neutral and ionised hy- 
drogen fractions like those presented in Fig. Al that we have 
performed more recently. 
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Figure Al. Test 1. Neutral and ionised hydrogen fractions obtained in simulations with the old (Pawlik &: Schaye 2008; first and third 
row) and new (second and fourth row) implementations of traphic. Shown are profiles of neutral and ionised fractions at times t = 30 
(left panel), 100 (middle panel) and 500 Myr (right panel), for simulations with (second and fourth row) and without (first and third row) 
resampling of the density field. The spatial resolution is fixed to Nqpu = 64^, iVngb = 32 and is indicated by the horizontal error bar in 
the upper left corner of each panel. The angular resolution is A^c = 128. The grey (light red) points show the neutral (ionised) hydrogen 
fraction for a randomly chosen subset of 10% of all particles. The solid black (red) curve shows the median neutral (ionised) hydrogen 
fraction in spherical bins and the error bars enclose 68.3% of the particles in each bin. The dashed black (red) curves show the exact 
solutions, obtained with our reference code ttId. The underestimate of the non-equilibrium neutral fraction exhibited in simulations 
with the old implementation of traphic is absent in the simulations that employ our new implementation, thanks to a new self-consistent 
manner of distributing photons absorbed by ViPs. The new implementation also reduces the scatter in the ionisation balance. 
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